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Preface 


The  Battlescale  Forecast  Model,  developed  at  the  U.S.  Army  Research 
Laboratory,  is  a  major  part  of  the  U.S.  Army  Integrated  Meteorological  System 
Block  II  software.  The  Battlescale  Forecast  Model  can  be  used  operationally 
world-wide  by  using  meteorological  data  obtained  through  the  Automated 
Weather  Distribution  System.  The  Integrated  Weather  Effects  Decision  Aids 
utilizes  the  Battlescale  Forecast  Model  output  file  to  calculate  the  values  of  their 
parameters.  The  Atmospheric  Sounding  Program  generates  products  such  as 
probability  of  thunderstorm  occurrence,  icing,  visibility,  and  turbulence. 

This  report  describes  the  theoretical  and  mathematical  principles  used  in  the 
Battlescale  Forecast  Model,  and  shows  some  examples  of  input/output  data. 
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Executive  Summary 
Introduction 

The  U.S.  Army  Research  Laboratory  (ARL),  Battlefield  Environment  Directorate 
(BED)  has  developed  the  Battlescale  Forecast  Model  (BFM)  for  operation  of 
short-range  (up  to  24-h)  forecasting  over  areas  extending  up  to  500  x  500  km. 
The  BFM  is  designed  so  it  can  be  applied  world-wide  as  long  as  meteorological 
data  is  available.  Depending  on  the  availability  of  the  input  data  from  the  U.S.  Air 
Force  Global  Spectral  Model  (GSM),  upper  air  sounding  data,  and  surface  data, 
different  combinations  of  input  data  are  used  to  initialize  the  BFM. 

Purpose 

This  report  describes  the  theoretical  and  mathematical  principles  used  in  the 
BFM.  The  GSM  forecast  data,  which  composes  a  major  part  of  the  input  data  to 
the  BFM,  is  compared  with  observed  upper  air  sounding  data.  The  influences  on 
the  BFM  forecast  field  of  surface  observation  data  assimilation  are  also  examined. 
Finally,  the  capability  to  forecast  stratiform  clouds  and  non-convective 
precipitation  by  the  BFM  are  demonstrated. 

Overview 

For  all  input  data  available,  the  BFM  is  initialized  by  the  composite  data  fields 
from  the  GSM  and  upper  air  sounding  data,  plus  surface  meteorological  data.  For 
a  minimal  requirement ,  the  BFM  can  be  initialized  by  a  single  sounding  profile 
within  the  model  domain  and  produce  forecast  field  for  6  h. 
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The  BFM  is  composed  of  elements  such  as: 

•  a  terrain  elevation  data  production  program; 

•  three  dimensional  data  analysis  program  for  model  initialization  and 
data  assimilation; 

•  a  prognostic  mesoscale  model,  the  Higher  Order  Turbulence  Model  for 
Atmospheric  Circulation  (HOTMAC); 

•  horizontal  and  vertical  displays  of  its  outputs; 

•  input/output  data  archiving  system;  and 

•  a  map  background  module  to  help  visualize  both  execution  and  display. 
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1.  Introduction 


The  U.S.  Army  Research  Laboratoiy  ,  Battlefield  Environment  Directorate  (ARL, 
BED)  has  developed  the  Battlescale  Forecast  Model  (BFM)  for  short-range 
weather  forecasting  over  areas  of  500  x  500  km  or  smaller.  The  BFM  can  also 
be  used  world-wide  by  using  meteorological  data  obtained  through  the  Automated 
Weather  Distribution  System  (AWDS).  It  is  a  major  part  of  the  U.S.  Army 
Integrated  Meteorological  System  (IMETS)  Block  II  software.  Meteorological 
data  includes  the  U.S.  Air  Force  Global  Spectral  Model  (GSM)  grided  forecast 
data,  in  addition  to  conventional  upper  air  sounding  and  surface  data. 

The  BFM  is  composed  of: 

•  A  terrain  elevation  data  production  program; 

•  Three-dimensional  data  analysis  programs  of  input  data  for  model 
initialization  and  data  assimilation; 

•  A  prognostic  mesoscale  model,  called  the  Higher  Order  Turbulence 
Model  of  Atmospheric  Circulation  (HOTMAC);  and 

•  Model  output  interpolation  programs  to  display  horizontally  (at  desired 
heights)  and  vertically  (at  model  levels)  for  any  desired  grid  location. 

The  BFM  has  been  used  successfully  in  exercises  such  as  the  Department  of 
Defense's  Operation  Desert  Capture  II  and  the  NATO's  Atlantic  Resolve  1994. 
Since  December  1995  it  has  been  employed  at  Heidelberg,  Germany,  by  the  617th 
Weather  Squadron  to  forecast  weather  over  Bosnia. 

The  objectives  of  this  report  are  to  describe  the  theoretical  and  mathematical 
principles  used  in  the  BFM  (mainly  elements  one  through  three),  and  to  show 
some  examples  of  input/output  data. 
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This  report  describes  the  method  of  obtaining  digital  terrain  elevation  data  and 
the  methods  of  three-dimensional  data  analysis  of  meteorological  data.  The 
prognostic  mesoscale  forecast  model,  HOTMAC,  is  summarized  from  the  original 
developer's  publications.  The  report  also  describes  the  methods  of  initialization 
and  data  assimilation.  Examples  of  the  BFM  input/output  are  given  to  show  the 
capabilities  of  the  model. 
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2.  Terrain  Elevation  Data  File  Production 


The  elevation  data  file  for  BFM  execution  is  produced  by  reading  the  Digital 
Terrain  Elevation  Data  (DTED)  CDROMs  from  the  Defense  Mapping  Agency 
(DMA).  DTED  is  a  uniform  matrix  of  digital  terrain  elevation  values.  Each 
matrix  point  represents  geographic  coordinates  with  a  corresponding  elevation 
value  above  mean  sea  level  (MSL).  Level  1  DTED  (DTEDl)  matrix  point  spacing 
is  3  arc  s  or  approximately  100  m.  [1] 

The  horizontal  grid  spacings  of  the  BFM  are  linear,  and  the  maximum  integration 
domain  has  a  horizontal  extent  of  500  km,  although  the  input  data  collection 
area  can  have  a  horizontal  extent  as  large  as  1600  km.  At  such  distances,  it  is 
a  good  assumption  that  the  earth  is  a  perfect  sphere  when  determining  the 
latitudes  and  longitudes  of  individual  grid  points.  However,  for  horizontal 
extents  beyond  about  1000  km,  this  assumption  loses  some  validity.  Still,  the 
errors  introduced  should  be  relatively  small. 

First,  we  calculate  the  latitudes  and  longitudes  of  the  BFM  grid  points  by  using 
the  following  equations.  Figure  1  shows  that  the  latitude  and  longitude  of  a  grid 
point  (i,  j)  is  calculated  as: 


GLON{ij)=LONQ  + 


Aii- 


^£'(T^)-cos[^-c;z^r(/)] 


180' 


180 


(1) 
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(2) 


GLAT(j)=LATQ+ 


(i  +1) 


180 


where 

Re  =  the  radius  of  the  earth 
A  =  the  unit  grid  space 

LONq  =  longitude  of  the  center  of  the  model  domain, 
LATq  =  latitude  of  the  center  of  the  model  domain. 


and 

(imax  jmax)  =  the  maximum  numbers  of  grid  points  for  the  (x,y)  directions, 
where 

(i,  j)  =  the  grid  point  counters  in  x-  and  y-  directions. 

The  values  of  i^^x  ^nd  are  fixed  in  the  model  to  values  of  161  x  161.  This 
means  that  for  a  horizontal  spacing  of  10  km,  terrain  data  must  be  generated  for 
an  area  of  1600  x  1600  km  for  the  three-dimensional  analysis  to  be  executed.  The 
southwest  comer  of  the  domain  is  represented  by  (1,  1),  and  the  northeast  comer 
is  represented  by  (i„3,^,  J. 
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LImax , 
,  J  m  a  X  J 
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(LON0, 

. LAT0 1 

n  .1) 

- ►  ^ 


Figure  1.  Schematic  diagram  showing  how  the 
latitude  and  longitude  of  model  grid  points  are 
calculated. 


Next,  the  terrain  height  at  (i,  j),  GZ  (i,  j),  is  interpolated  from  the  data  of  the  four 
surrounding  matrix  points  that  are  read  from  the  DTED  CDROM.  Referring  to 
figure  2  let’s  define: 

GLON{ij)-LON^ 

^  (3) 

LON-LON^ 

GLAT{i)-LAT^ 

u= -  (4) 

LAT^-LAT^  ^  ’ 
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Here, 


t  —  the  dimensionless  longitudinal  distance  to  (i,  j)  from  the  western  longitude 
of  the  DIED  data  that  surround  the  grid  point  (i,  j). 

Similarly, 

u  =  the  dimensionless  latitudinal  distance  to  (i,  j)  from  the  southern  latitude 
of  surrounding  DTED  data. 


2 

^x  +  l  ,  y  +  1  1 

GZI i .j 
L.ONw  ° 

L 

1  fl  T  _ _ 

) 

DNe 

2 

4 

r  t  1 5 

[(x.yl 

_ 

nx  +  1  ,y) 

1 - r 

Figure  2.  Schematic  diagram  showing  how 
the  elevation  values  are  obtained  from  the 
DTED  CD-ROM  Level  1  data. 

The  terrain  elevation  at  the  grid  point  (i,  j)  is  determined  by 

GZ(i,  j)  =  ( 1  -t)-(  1  -  u)-  Z(x,  y)  + 1-  ( 1  -  u)-  Z(x  +  1 ,  y)  (5) 

+  (1  - 1)-  u-  Z(x,  y  +  1)  +  f  u-  Z(x  +  1,  y  +  1) 
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where 


GZ  =  terrain  elevation  of  the  BFM  grid  point, 
Z  =  a  DTED  elevation  data  point. 


The  above  process  is  repeated  for  eveiy  grid  point  of  the  BFM. 


Figure  3  is  an  example  of  elevation  data  obtained  for  the  southwestern  United 
States  centered  at  the  National  Training  Center,  Fort  Irwin,  CA.  Data  is  produced 
for  the  area  of  1,600  x  1,600  km  with  unit  grid  spacing  of  10  km.  The  internal 
square  area  is  the  BFM  domain,  covering  500  x  500  km.  Both  grid  squares  share 
the  same  center  coordinates. 


i  ctmus 


Figure  3.  Terrain  elevation  data  produced 
over  southern  California,  centered  at  the 
National  Training  Center,  Fort  Irwin,  CA. 
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3.  Three-Dimensional  Analysis  of  Input  Data 


The  BFM  forecast  calculations  are  made  by  using  the  GSM  grided  forecast 
fields,  upper  air  radiosonde  observations,  and  surface  sensor  observations  as 
initial  input  and  time-dependent  lateral  boundaiy  condition  data.  This  data  is 
regularly  communicated  to  the  IMETS  workstation  computer  through  the  AWDS 
from  the  Air  Force  Global  Weather  Center  (AFGWC),  or  manually  input. 
Objectively  analyzed  three-dimensional  data  fields  for  model  initialization  can  be 
provided  by  three  different  data  sources: 

•  GSM  data, 

•  Upper  air  data,  and 

•  Combination  of  GSM  and  upper  air  data. 

Each  of  the  above  can  be  considered  a  separate  analysis  “mode.”  Depending  on 
the  availability  and  the  quality  of  data,  any  of  the  above  modes  can  be  selected  for 
initialization  analysis.  Surface  sensor  data  is  not  used  to  make  three-dimensional 
data  fields;  it  is  used  to  initialize  the  BFM  prognostic  calculations  in  a  different 
way. 

Figure  4  shows  the  distribution  of  the  different  input  data.  "G"  represents  the 
GSM  data  locations,  "U"  the  upper  air  radiosondes,  and  "S"  the  surface  data  sites. 
The  center  of  the  model  domain  was  chosen  to  be  at  18.657  E  and  44.470  N.  The 
entire  domain  covers  a  1600  x  1600  km  area  (161  x  161  grid  points  with  a  unit 
grid  distance  of  10  km),  and  the  inner  square  covers  a  sub-area  of  500  x  500  km 
(51x51  grid  points).  The  large  square  covers  a  1,600  x  1,600  km  area.  The  inner 
square  of  the  BFM  domain  covers  a  500  x  500  km  area.  The  area  is  centered  at 
Tuzla,  Bosnia,  located  at  44.470  N  and  18.657  E.  “G”  is  located  at  the  GSM  data 
points,  “U”  at  the  upper  air  sounding  sites,  and  “S”  at  the  surface  observation 
sites. 
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Figure  4.  Distribution  of  input  data  over 
the  BFM  domain. 

The  GSM  data  and/or  upper  air  radiosonde  data  available  in  the  large  square  is 
used  to  produce  objectively  analyzed  three-dimensional  meteorological  data  fields 
at  the  horizontal  resolution  required  by  the  model.  Then,  the  data  fields  within  the 
inner  square  sub-region  are  used  to  initialize  the  BFM  forecast  model. 

In  all  the  above  three  methods,  the  analysis  was  performed  over  a  large  domain 
of  161  X  161  grid  by  using  z*  defined  by  equation  (6),  in  which  is  the 
maximum  height  in  the  161  x  161  grid  domain.  A  forecast  calculation  is 
typically  done  over  a  51  x  51  grid  domain  located  at  the  central  part  of  the  large 
domain  where  the  maximum  height  is  z^^  .  In  the  case  where  Zg^^,,  is  smaller  than 
Zgmax2»  all  the  parameters  are  adjusted  linearly  to  the  new  z*  levels  defined  by 
replacing  Zgmax-  The  three-dimensional  fields  covering  the  51  x  5 1  x  16 

subgrid  are  used  to  initialize  the  BFM. 
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3.1  GSM  Data 


The  GSM  forecast  calculations  are  made  regularly  at  the  base  times  of  00  and 
12  Greenwich  Meridian  Time  (GMT).  Forecast  results  become  available  to  the 
IMETS  approximately  7  h  after  the  base  times.  For  instance,  if  the  BFM  forecast 
calculation  is  made  at  a  base  time  between  00  and  12  GMT,  the  GSM  forecast 
data  based  at  00  GMT  is  not  available  until  about  07  GMT;  the  calculations  are 
based  at  12  GMT  of  the  previous  day  are  available.  To  avoid  the  data-absent 
period,  the  12,  24,  36,  and  48-h  forecast  data  generated  at  the  base  time  of 
12  GMT  of  the  previous  day  are  used.  When  the  BFM  forecast  base  time  is 
between  the  12  and  24-h  GSM  forecast  time,  the  initial  BFM  field  is  constructed 
by  linearly  interpolating  between  the  12  and  24-h  GSM  forecast  fields  (in  addition 
to  combining  and  compositing  available  upper  air  radiosonde  data).  Similarly,  the 
GSM  derived  forecast  fields  used  by  the  BFM  as  time  dependent  lateral  boundary 
values  (out  to  24  h),  can  be  obtained  by  linear  interpolations  of  the  appropriate 
GSM  fields:  24,36,  and  48  h. 

Due  to  current  design  constraints,  although  a  new  GSM  base  time  forecast  arrives 
to  the  IMETS  about  7  h  later,  this  data  is  not  accessed  until  exactly  12  h  after  the 
new  base  time.  Therefore,  for  a  BFM  initialization  1 1  h  after  either  the  12  GMT 
or  00  GMT  GSM  base  times  (23  GMT  or  1 1  GMT,  respectively),  GSM  23-h 
forecast  data  is  used  to  generate  the  BFM  initialization  fields  (along  with 
available  upper  air  and  surface  data).  Additionally,  GSM  47-h  forecast  data 
would  be  used  as  the  BFM  24-h  lateral  boundary  condition  data.  A  new  design 
is  being  considered  that  would  make  the  new  GSM  base  time  forecasts  available 
to  the  BFM  immediately  upon  arrival  to  the  IMETS. 

The  following  GSM  products  on  the  whole-mesh  grid  (unit  grid  space  is  381  km 
at  60°  N  or  S  and  a  function  of  latitude  [2]  )  on  six  constant  pressure  levels  (200, 
300,  500,  700,  850,  and  1000  mb)  are  used: 
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=  Geopotential  height, 

T  =  Temperature, 

(T  -  Tj)  =  Dew  point  depression, 

(U,  V)  =  Horizontal  wind  vector  components. 

A  fictitious  1070  mb  level  is  created,  with  data  extrapolated  down  from  the 
1000  mb  level  to  ensure  that  data  exists  down  to  the  minimum  terrain  elevation 
of  the  161  X  161  grid.  Using  the  above  parameters,  three-dimensional  fields  on 
(x,  y,  z*)  coordinates  of  the  following  parameters  are  produced: 

©v  =  Virtual  potential  temperature  (K), 

Qy  =  Water  vapor  mixing  ratio  (g/kg), 

U  =  East-west  component  of  wind, 

V  =  North-south  component  of  wind. 

Pgr  =  Surface  pressure  distribution, 

where 

z*  =  vertical  coordinate  used  in  HOTMAC,  and  defined  as: 


_z-z„ 
=H- — ^ 


H-z. 


(6) 


where 

z  =  the  Cartesian  vertical  coordinate, 

^  =  the  ground  elevation, 

H  =  the  material  surface  top  of  the  model  in  the  z*  coordinate. 
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and 


H  =  the  corresponding  height  in  the  z  coordinate  defined  by  H=H  + 


where 

Zgmax  ^  the  maximum  value  of  the  terrain  elevation  in  the  BFM  domain. 

Currently,  the  vertical  extent  of  the  model  is  7000  m  above  the  highest  grid  point 
of  the  BFM  domain.  Table  1  show  the  variables  as  calculated  at  the  following 
16  vertical  levels: 


Table  1.  Z*  heights  for  vertical  levels  of  BFM 


Level 

(m) 

z*  height 

1 

0.0 

2 

2.0 

3 

6.0 

4 

10.0 

5 

14.0 

6 

32.3 

7 

151.0 

8 

384.5 

9 

732.6 

10 

1195.4 

11 

1773.0 

12 

2465.3 

13 

3272.2 

14 

4193.9 

15 

5230.3 

16 

6381.4 
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Figure  5.  Model  vertical  level  distribution 
over  complex  terrain. 


Figure  5  shows  the  model  vertical  level  distribution  over  complex  terrain, 
showing  that  above  a  higher  location,  vertical  layers  are  more  densely  packed. 
Layers  near  the  surface  are  not  shown. 

Here, 

z  =  the  model  level  number. 

The  calculation  of  z*  is  z*(k)=  c,(k-0.5)^  +  CjCk-O.S)  +  C3, 
where 

k  =  the  level  number, 
and 

c„  C2,  and  C3  =  constants  determined  within  the  model. 
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Figure  5  is  a  schematic  presentation  of  how  the  BFM  vertical  layers  are 
distributed  over  complex  terrain.  The  layers  are  more  densely  packed  over  a  high 
peak  than  over  a  low  elevation  point.  More  layers  exist  near  the  surface  of  the 
earth. 

The  following  procedures  are  taken  to  analyze  the  data: 

1.  Horizontal  wind  components  U' and  V  (on  the  whole-mesh  reference  mesh 
of  the  GSM)  are  converted  to  U  and  V,  where  U  and  V  are  (respectively)  the 
east-west  and  the  north-south  components  of  horizontal  wind  vectors. 

By  defining 


A' =  27c  -  [(7r/180)-80  +  A] 


(7) 


where 

A  =  the  longitude  of  the  GSM  grid  (positive  for  the  eastern  hemisphere). 
U  and  V  are  calculated  by: 


U  =  U'-cos  A’  -  V'-sin  A' 

(8) 

V  =  U'-sin  A’  +  V'-cosA' . 

(9) 

2.  Water  vapor  mixing  ratio,  Q^,  at  the  pressure  level,  P,  is  calculated  from 
dewpoint  temperature: 


Q  =0.622- 
P 


(10) 
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(11) 


e=6.11exp[ 


MllT, 

- -] 

Tj+2313 


where 

Tj  =  the  dew  point  temperature  in  °C  calculated  as  =  T  -  A, 
A  =  the  dew  point  depression, 
e  =  the  saturation  vapor  pressure  at  the  temperature  Tj 

and 


P=  the  pressure. 

Virtual  temperature,  T^,  is  calculated  as: 


T,  =  T(1+O.610v). 


(12) 


3.  The  parameters,  U,  V,  T,  and  O  on  pressure  levels  are  horizontally 
interpolated  to  the  BFM  horizontal  grid  points  by  the  following  method.  A  first 
guess  value  of  a  parameter  Y  at  a  grid  point  (i,  j)  is  calculated  as: 


Y(y> 


4k 


Ea^cxpC 


2 

Ak^ 


(13) 
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Here 


Y  =  any  parameter, 

exponential  function  =  the  Barnes's  weighting  function,  [  3  ] 

rjjN  =  the  normalized  distance  between  a  grid  point  (i,  j)  and  the  Nth  GSM 

point, 

and 

k  =  an  empirical  parameter  to  determine  the  shape  of  the  weighting  function. 

By  using  the  first-guess  values  of  the  four  grid  points  surroxmding  the  Nth  GSM 
point,  an  interpolated  value  at  that  point,  ,  is  bilinearly  calculated  as: 

tj  =  T(x,  y)  +  (x*  -  x)-[Y(x  +  1 ,  y)  -T(x,  y)]  (14) 

t2  =  Yfx,  y+1)  +  (x'  -  x)'[Y(x  +  1,  y  +  1)  -  Y(x,  y  +  1)  (15) 

Y'n=  ti  +  (y'-y)’[i2 -til-  0^ 


Here, 

(x,y)  =  the  southwest  grid  of  four  grid  points  surrounding  the  Nth  GSM 
point  located  at  (x',y'). 

In  figure  6, 


"  *  "  =  the  locations  of  GSM  data  points. 
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Figure  6.  Schematic  diagram  showing  the 
relationships  between  the  BFM  grid  points 
and  GSM  data  points. 


The  difference  between  the  GSM  value  Y  ^  and  Y'  ^  ,  A ,,,  =  Y  ^  -  Y'  n, 
distributed  to  the  BFM  grid  points  as: 


EvVxp(-^) 


A(y)=- 


4y 


E.exp(-^) 

4y 


where 


Y  =  an  empirical  weight  reduction  factor  (0.2). 


is  now 


(17) 


A  final  interpolated  value  of  Y  at  (i,  j),  Y/i,  j)  is  calculated  as: 


TAj)  =  W,j)  +  A(i,j).  (18) 

4  The  parameter  Y  on  pressure  surfaces  are  linearly  interpolated  to  z*  level  of 
the  BFM  as: 


Y  -Y 


(19) 


where 


=  the  Cartesian  height  above  sea  level  of  z*,  calculated  from  equation  (6) 
as: 


z  =z  +z 

si  g 


(H+z  -z ) 

gmax  ^gf 

H 


(20) 


5.  The  surface  pressure  distribution  over  a  complex  terrain  is  determined  as 
follows.  If  the  terrain  height  Zg  is  between  the  pressure  surfaces  P,  and  Pj  (Pj  > 
Pj )  where  the  geopotential  heights  and  the  temperature  are,  respectively,  (j) ,  and 
(4>i  )  ^“i<l  T  ,  and  T  the  surface  pressure  is  calculated  using  the 

hypsometric  equation  as: 

I)Zg-(l)i  <  (j>2-Zg 


where 


and 


ii)  z,-<|),>(|>l,-z. 


where 


and 


T.-T,  _ 
4)2-q>i 


(22) 


Zg^-01>$2-Z^ 


(23) 


^^^=^2exp[-^(c|)2-0] 


RJ 


(24) 


r,-?,  _ 

r=7-,-^^(z-z,) 

Z2  Zj 


(25) 


^=|(^^+4)2) 


(26) 
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3.2  Upper  Air  Sounding  Data 


Prior  to  the  three-dimensional  objective  analysis,  a  quality  control  of  upper  air 

sounding  data  is  made.  The  following  criteria  are  applied  to  each  sounding  data: 

.  Flag  the  geopotential  height  values  if  the  values  differ  by  1 5  percent  or  more 
from  the  standard  atmosphere  for  pressure  levels  between  500  and  880  mb, 
and  differ  by  20  percent  or  more  from  the  standard  atmosphere  for  pressure 
levels  below  880  mb.  This  check  is  not  done  at  the  first  two  levels  above  the 
surface. 

•  Flag  the  height  of  the  second  level  if  an  extremely  small  pressure  gradient 
(<0.02  mb/m)  near  the  surface  is  found. 

.  Eliminate  all  layers  with  extreme  temperature  inversions  (AT/Az  > 
4°C/  100  m),  except  near  the  surface.  If  such  an  inversion  exists,  a  lower 
level  is  flagged.  Near  the  surface,  an  inversion  as  extreme  as  400°  C/km 
is  allowed. 

.  Eliminate  extremely-superadiabatic  layers  (AT/Az:^-22  °C/100  m)  by  flagging 
the  lower  level,  except  near  the  surface  where  lapse  rates  less  than  the 
autoconvective  lapse  rate  (-  g/R^  =  -3.42°C/100  m)  are  allowed.  Those 
critical  values  in  the  above  are  arbitrary,  but  seem  to  be  not  too  unrealistic. 
More  rigid  criteria  seem  to  often  eliminate  valid  temperature  data  from 
observed  data  which  could  be  of  importance  to  modeling  the  boundary  layer. 

•  Eliminate  vertical  wind  speed  shears  exceeding  150  kn/km  by  flagging  both 
levels. 
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.  Correct  duplicate  pressure  levels  or  heights/pressure  out  of  order. 

.  Edit  the  data  manually  from  the  GUI  window. 

After  the  above  steps,  check  the  height  fields  by  using  the  hypsometric  equation. 

After  the  above  quality  control,  take  the  following  procedures  to  produce  three- 

dimensional  data  fields  from  the  upper  air  sounding  data: 

.  At  each  sounding  location,  parameters  such  as  horizontal  wind  vector 
components,  potential  temperature,  and  water  vapor  mixing  ratio  are  vertically 
interpolated  to  20  different  height  levels  (Zj  ASL)  by  using  a  linear 
interpolation  method. 

•  At  each  Z;  level,  horizontal  interpolation  to  the  BFM  grid  is  performed  by 
using  the  weighting  factor,  1/r^  as: 


^  2 

V  — 

2 

^ijJ^ 


(27) 


The  weighting  factor  1/r^  is  used  here.  According  to  our  experience,  for 
a  small  number  of  observations,  the  Barnes  method  does  not  seem  to 
produce  reasonable  fields. 

.  Linear  vertical  interpolation  from  z  ^  to  z  given  by  equation  (20)  is 
performed  for  each  parameter. 
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•  The  pressure  values  at  H  „jt  m  height  ASL,  P  Hcrit,N  the  sounding  site  N  is 
calculated  for  each  sounding  site.  The  horizontal  distribution  at  H  jn,  m  level, 
PHcrit(IJ)  is  calculated  by  using  the  weighting  factor  1/r^ .  The  level  of 
m  is  chosen  so  that  almost  all  sounding  data  have  observed  values.  H  can  be 
either  100,  500, 1000,  1500, 2500,  or  3000  m  ASL,  and  is  determined  from 
the  available  radiosonde  site  elevations. 

.  The  difference  of  pressure  between  P  Hen,,  n  and  the  nearest  grid  value  of  P  henp 
A  N,  is  redistributed  horizontally  to  grid  points  by  using  the  same  weighting 
factor,  and  the  new  pressure  values  at  H„jt  m  level,  P^gri,  are  calculated  as: 

Y' 

2 

^  (28) 

y  — 

2 

r 

^  ij^ 

.  The  surface  pressure  Pgn,d  is  calculated  from  PHcnthy  using  the 
hypsometric  equation: 

-ZgiU))]  ^29) 

RJj 


where 


T  =  the  average  absolute  temperature  between  Zg  and  m  level,  and 
inside  the  exponential, 
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+  sign  -  Zg  below  meters  ASL 
-sign  =  Zg  above  H^ritn^eters  ASL. 


.3  Compositing  the  GSM  and  Upper  Air  Sounding  Data 

When  both  the  GSM  data  and  upper  air  sounding  data  are  available,  three- 
dimensional  fields  of  the  parameters  composited  fi-om  both  data  sources  are 
usually  used  to  initialize  the  BFM.  The  following  procedures  are  taken  to 
composite  the  two  different  data: 


•  Three-dimensional  fields  of  the  parameters  (virtual  potential  temperature, 
dew-point  temperature,  horizontal  wind  components,  u  and  v)  are  obtained 
fi'om  the  GSM  data  on  20  vertical  layers,  Zj, 

.  Upper  air  sounding  data  is  also  interpolated  to  the  20  vertical  layers,  z  j. 


.  Three-dimensional  data  fields  created  fi-om  the  GSM  data  are  now  used  as 
background  fields,  T Q(i,  j),  and  composited  with  the  upper  air  sounding  dat^ 
Yu.  Let's  define  the  mean  difference,  d  as: 


E 

d=^ 


N 


(30) 


where 

N  ^ke  interpolated  value  to  the  Nth  upper-air  location  by  using  the 
four-surrounding  grid  points,  and  Yu,n  Nth  upper-air  data. 


A  bilinear  interpolation  method  described  in  equations  (14),  (15)  and  (16)  is 
used  to  calculate  Yq  n. 

Byreplacing  YG(i,j)  by 

n(iJ)  =  ’i'G(iJ)+d  (31) 

the  mean  error  is  removed,  but  the  values  of  n  interpolated  from  Y*G(i,  j)  for 
the  Nth  upper  air  location  are  not  generally  equal  to  Yuj^. 

The  difference 

(i*  =W  _iir* 

“  N  ^  U,N  ^  G,N 

is  not  zero,  but  the  mean  d*  is: 

7* 

N 


N 

V=1 


(32) 


(33) 


The  next  step  is  to  distribute  d*  to  the  entire  field  using  a  weighting  function  of 
1/r^as: 


E 

E 


N=\ 


2 

1 


iJ,N 


(34) 
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The  above  procedures  are  applied  to  all  vertical  levels,  Zj,  for  all  the  parameters, 
and  vertical  linear  interpolation  from  Zjto  z*  are  performed.  Linear  interpolation 
from  Zj  levels  to  z  5,  are  performed  for  all  the  parameters,  and  dew-point 
temperature  field  is  converted  to  water  vapor  mixing  ratio  field  by  using  the 
formula  given  in  equation  (11). 

For  surface  pressure  calculations,  the  observed  surface  pressure  values  are 
corrected  to  the  model  elevation  by  using  the  hyposometric  equations  because 
the  radiosonde  site  elevations  and  the  model  elevations  for  the  corresponding 
locations  are  generally  different.  The  corrected  observed  values  are  composited 
by  using  the  above  procedures. 
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4.  Forecast  Model 


For  forecast  calculations,  the  BFM  uses  a  mesoscale  model  (HOTMAC) 
developed  by  Yamada.  [4]  This  model  is  based  on  a  set  of  second-moment 
turbulence  equations  closed  by  assuming  certain  relationships  between  unknown 
higher-order  turbulence  moments,  and  the  known  lower-order  moments.  The 
model  output  variables  are  wind,  potential  temperature,  the  mixing  ratios  of 
water  vapor  and  cloud  water,  turbulence  second  moments,  a  turbulence  length 
scale,  and  turbulence  transport  coefficients  (eddy  viscosity  and  eddy  diffusivity). 
The  model  is  time  dependent  and  three-dimensional  in  space. 

HOTMAC  can  be  used  under  quite  general  conditions  of  flow  and  thermal 
stratification.  Methods  for  turbulence  parameterization  are  more  advanced  than 
those  in  simple  eddy  viscosity  models.  This  model,  combined  with  a  statistical 
cloud  model,  can  simulate  interaction  between  water  phase  changes  and  basic 
dynamic  variables.  Effects  of  short  and  long-wave  solar  radiation,  and 
topography  are  also  included  in  the  model.  Surface  temperature  is  computed 
fi’om  a  heat  conduction  equation  for  the  soil,  and  a  heat  energy  balance  equation 
at  the  surface. 

The  model  assumes  hydrostatic  equilibium  and  uses  the  Boussinesq 
approximation.  Under  the  Boussinesq  approximation,  it  is  assumed  that  the 
modeled  fluid  is  incompressible  to  the  extent  that  thermal  expansion  produces 
a  buoyancy.  Buoyancy  forces  are  retained  in  a  hydrostatic  basic  state  with  respect 
to  pressure  and  density  via  the  inclusion  of  small  pressure  and  density  deviations. 
This  assumption  holds  as  long  as  the  density  deviation  from  the  mean  state  is 
negligible.  The  density  variations  are  only  considered  when  they  are  closely 
coupled  to  gravity.  Therefore,  in  theory,  the  model  applications  are  limited  to 
flows  where  the  local  acceleration  and  advection  terms  in  the  equation  of  vertical 
motion  are  much  smaller  than  the  acceleration  due  to  gravity  (hydrostatic 
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equilibrium)  and  temperature  variations  in  the  horizontal  directions  are  not  too 
large. 

The  content  of  this  chapter  is  mainly  summarized  from  publications  by  Yamada. 
[3]j  [4]  ,  [5]  Some  modifications  are  made  to  the  original  source  code,  and 
noted  in  the  text. 


4.1  Governing  Equations 


A  terrain  following  vertical  coordinate  system  is  used  to  increase  the  accuracy 
in  the  treatment  of  surface  boundary  conditions  equation  (6).  Following  a 
coordinate  transformation  [5],  the  governing  equations  are  given  by: 


DU  H-r  *  <Q>  dz 

Dt  ^  0.  dx 


dx  "dx  dy  "^dy  H-z 


dU^  ,  H  a  ,  — 


dz* 


(-UW) 


(35) 


dx  "^dx  H-z 


^  (-vw) 


dz  * 


(36) 


where 


dU ^dV  ^  dW* 
dx  dy  dz  ’ 


1  dz  dz^ 

H-z^  dx  dy 


(37) 


W*=- 


H 


H-z 


z  *  -H  dz^  dz 
H-z  dx  dy 


(38) 
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and 


Dt  dt  dx  dy  Sz  * 


(39) 


In  the  above  expressions,  <0v  >  is  the  base  state  virtual  potential  temperature, 
which  we  are  defining  as  the  initial  virtual  potential  temperature  field  obtained 
by  the  three-dimensional  objective  analysis  of  input  data.  The  perturbation 
potential  temperature  field  is  defined  as  the  change  fi'om  this  base  state  over 
time.  Due  to  this  approach,  in  cases  with  very  strong  fi’ontal  zones  and  very 
large  temperature  gradients,  the  Boussinesq  approximation  may  lose  some  of  its 
validity.  The  second  terms  on  the  right-hand  side  of  equations  (33)  and  (34)  and 
the  fourth  term  of  equation  (35)  express  the  effects  of  ground  slope. 


In  the  original  code  of  the  HOTMAC,  the  geostrophic  winds  U  g  and  Vg  are 
computed  by: 


H 


<©  >  dy 


Hdyh- 


<©  > 


*  *  IJ  Jz- 


<QIH)>  H  <Q>dx 


Hdxh^  <©> 


(40) 

(41) 


where 


60  =0.-  <0v>, 


and  the  abbreviated  symbols 
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Ug(H)=Ug(x,  y,  H,  t) , 


and 


Vg(H)=Vg(x,  y,  H,  t) 


are  used. 

In  deriving  equations  (38)  and  (39),  Yamada  makes  several  assumptions,  the 
most  substantial  being  that  <0^>  is  independent  of  height  z.  [5]  Based  on  our 
definition  of  <0^>,  this  is  only  a  valid  asssumption  in  the  lowest  1  to  2  km  AGL. 
It  is  invalid  for  the  BFM  since  its  model  top  is  at  7000  m  AGL  or  greater. 
However,  by  eliminating  this  assumption,  we  have  detected  that  equations  (38) 
and  (39)  produce  unrealistically  extreme  magnitudes  of  geostrophic  wind  speeds, 
particularly  over  complex  terrain.  This  is  due  to  the  problem  of  trying  to  resolve 
the  difference  between  two  very  large  terms  of  opposite  sign,  needed  to  compute 
the  pressure  gradient  in  the  z*  coordinate  system.  To  avoid  this,  the  following 
approach  has  been  taken: 

•  The  values  of  <0^  (H)>,  Ug(H),  and  Vg  (H)  in  equations  (38)  and  (39)  are 
replaced  by  the  values  at  the  grid  point  vertical  level  nearest  to 
1500  mAGL, 


•  Below  level  z.^,  solve  for  Ug  and  Vg  using  equations  (38)  and  (39)  and  the 
substitutions  mentioned  above  in  1 . 

•  At,  and  above  level  Ug  and  V  g  are  set  equal  to  the  smoothed  model 
generated  grid  point  wind  value,  using  the  four  surrounding  grid  points  and 
a  smoothing  factor  of  0.5.  At  the  lateral  boundaries  the  same  scheme  is 
applied,  except  that  the  three  surrounding  grid  points  are  used. 
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The  above  process  eliminates  the  false  generation  of  extreme  Ug  and  Vg  wind 
components ,  and  allows  the  large-scale  geostrophic  components  to  stay  in  line 
with  the  large-scale  GSM  guidance.  At  the  same  time,  it  allows  the  geostrophic 
components  in  the  lower  boundary  layer  to  be  determined  by  the  pressure 
gradient  force  (using  thermal  wind  principles  in  the  layer  z  *  to  and  the  model 
perturbation  virtual  potential  temperature  field)  generated  by  the  BFM  surface 
physics.  The  geostrophic  wind  calculating  scheme  may  be  modified  in  the  future. 


A  turbulence  kinetic  energy  equation  is  given  by: 


H 


(^)] 


d 


H-z  dz 


dz*  2 


H  ,—  dU  —dV.  „ 

{UW - +VW - )+PgWD  -  ^ 


(42) 


H-z^  dz*  dz 
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a  turbulence  length  scale  /  is  obtained  from: 


H-z^  dz 


dz* 


H  .  —dU  —dV.  o  -5-, 

[77 ( -  ww - -MW - ) +P^6  J 

H-z  dz*  dz* 


(43) 


where 

=  twice  the  turbulence  kinetic  energy  defined  as: 

2  2  2  9 


(44) 


where 

w6y  =  turbulence  heat  flux, 

0v  =  the  fluctuation  part  of  virtual  potential  temperature. 


and 


(F,_F2,  Sq,  S,  B)=(L8,  1.33,  0.2,  0.2,  16.6.). 
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These  empirical  constants  were  determined  from  laboratory  experiments  [7]  for 
neutral  boundary  layer  flow,  and  are  applied,  in  the  BFM,  to  other  atmospheric 
stability  cases.  The  internal  heat  energy  equation  is  written  as: 


Db&  ^  fj  50  ij 

— — I]  +  A[/i:  _ji] 

Dt  dx  dx  dy  ^  dy 


dz* 


1 

pCp  dz  * 


a00„ 

-w - ^ 

dz 


(45) 


The  long-wave  radiation  flux  R^/pCp  is  computed  according  to  Sasamori.  [8]  A 
conservation  equation  for  the  mixing  ratio  of  total  water  (Q  w=Qv  +  Qc)  is  given 
by 


Dt  dx  ^  dx  dy  ^  dy 


(46) 


where 

Qy  =  the  mixing  ratio  of  water  vapor 
and 


Qc  =  the  mixing  ratio  of  cloud  water. 


The  turbulent  fluxes  in  equations  (33),  (34),  (41),  (42),  (43),  and  (44)  are 
obtained  from  simplified,  second-moment,  turbulence-closure  equations,  [9] 


41 


{uw,vw)=-lqSJ_^~] 
oz  dz 


(47) 


(wQ,wq^)=-alqSJ^^,^] 

OZ  OZ 


(48) 


where 

Sm  and  a  =  functions  of  the  flux  Richardson  number, 
a  is  the  reciprocal  of  the  turbulent  Prandtl  number  and  given  by  K  h/K^ 
where 

Kjj  =  an  eddy  diffusivity  coefficient 
and 

Km  =  an  eddy  viscosity  coefficient. 


The  final  expressions  for  and  a  are  given  in  Yamada  [9]  as: 
forRif<0.16 

5  =1  n^(Q-1912-ig9(0.2341-ig/^) 

^  ■  (I -RiyXO. 2231 -Rij) 


(49) 
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(50) 


and  for  Rif  >  0.16, 

Sm  =  0.085 


0.2231-/?/, 

a=1.318 - ^ 

0.2341  -Rij. 


Rif  <  0.16 


(51) 


=  1.12  Rif  >0.16 

The  relationship  between  the  gradient  Richardson  number 

_  gdQldz 
^  QidUIdzf 


(52) 


and  the  flux  Richardson  number 


gwQ 

Qii^dUIdz+^dVIdz] 


(53) 


is  given  by: 

/?z}.=0.6588(/e/^+0. 1 116-[Rf  -0.322 1/?/+0.03 1 56]’'^) 


for  R,g<Ri,  =  Rif,  for  R,g  >Ri, 


43 


where 


~  ^  critical  flux  Richardson  number 


and 


Ric  (-0.195)  -  a  critical  gradient  Richardson  number. 

4.2  Surface  Boundary  Conditions 

Surface  boundary  conditions  for  equations  (35),  (36),  and  (42)  through  (46)  are 
constructed  from  the  empirical  formulas  by  Dyer  and  Hicks  for  the 
nondimensional  wind  and  temperature  profiles  as  follows:  [10] 

(55) 


e(z)-9(0) 


=-^[ln(^)-T,(0] 


(56) 


Q. 


■=f[ln(-^)-Y/0 


-Ov 


(57) 


(58) 


q\z)^kzq\z) 


(59) 
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In  the  above  equations,  V(z),  ©(z)  and  Q  v(z)  are  abbreviations  for  V(x,y,z,t),  0 
(x,  y,  z,  t),  and  Qy(x,  y,  z,  t),  respectively,  and  V(z)  is  the  horizontal  wind  speed. 
Terms  u„  T.  and  Q,  defined  as  u  .=(T/p)'^,  T  ,=H/pc  pU »,  and  Q  ♦=  E/pu . ,  are 
friction  velocity  and  scales  for  temperature  and  water  vapor,  respectively,  t,  H,  and 
E  are  surface  stress,  total  sensible  heat,  and  rate  of  evaporation,  respectively.  Total 
sensible  heat  is  defined  as: 

=w0y=(l+O.61g^)w0+O.61w0  (60) 


and  the  Monin-Obukhov  length,  L,  is: 
L=-ulll^gH 


(61) 


The  parameter 

C  =  a  nondimensional  height  z/L; 
k  =  the  von  Karman  constant; 

Zq,  Zot,  and  Zqv  =  roughness  lengths  for  wind,  temperature,  and  water  vapor; 
and 

Pj  and  S(.  =  turbulence  Prandtl  and  Schmidt  numbers,  respectively. 

Terms 

T„,T,,andYJll]  =  correction  terms  for  the  atmospheric  stability. 

Their  functional  forms  are  given  by: 
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(62) 


d 


(63) 


(64) 


where 

4>m5  4>h’  ^  nondimensional  wind,  temperature,  and  water  vapor  gradients. 

The  following  formulations  [10]  are  used  for  (j)„,  4)^,  and  (j)^  under  unstable 
conditions: 


<|)„(C)Ka./fc)(a£//az)>(l-I5C)-'« 

<l)/0=(7'yfe)oe/aj)=(i  -i50‘“ 

<t>,(0  =(e./fexae/az) =(  1  - 1 50''° 


(65) 

(66) 
(67) 


For  stable  conditions. 


(l)„(C)=c|),(0=(t>,(0-i-^5C 


(68) 
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The  above  formulas  are  valid  only  for  horizontally  homogeneous  surfaces. 
However,  it  is  assumed  that  the  same  relations  are  fair  approximations  over 
nonhomogeneous  terrain,  provided  that  the  formulas  are  applied  sufficiently  close 
to  the  surface.  Vegetation  plays  an  active  part  in  the  apportionment  of  available 
heat  energy  between  convective  (sensible  and  latent)  and  conductive  (into  soil) 
components.  HOTMAC  has  a  parameterization  scheme  of  tall  tree  effects,  but  in 
the  current  BFM  it  is  assumed  that  an  entire  model  domain  is  covered  with  bare 
soil.  In  the  future,  the  parameterization  ofvegetation  cover  will  be  included  in  the 
model. 


In  the  model,  the  distinctions  between  land  and  water  are  made  in  surface 
roughness  parameters  and  surface  temperature.  Over  land,  Zq  and  Zq,  (=Zov)  are, 
respectively,  0.1  and  0.0135  m.  Over  water,  Zo=0.016u.Vg,  and 
Zot  =  Zoy  =  0.000022/ku.  Surface  temperature  over  land  is  calculated  from  the 
surface  energy  balance.  Temperature  over  water  is  assumed  to  stay  constant 
during  a  forecast  period  at  the  values  determined  from  the  climatological  monthly 
average  temperatures. 

4.3  Surface  Energy  Balance  and  Surface  Temperature 


Use  of  the  similarity  formulas  requires  knowledge  of  surface  temperatures.  The 
temperature  Tj  in  the  soil  layer  is  obtained  by  solving  the  heat  conduction 
equation: 


a 

dz/  "az 


(69) 


where 
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Zs= positive  downward, 


and  soil  heat  conductivity  can  be  a  function  of  soil  moisture  content  and 
emissivity  e.  Appropriate  boundary  conditions  for  solution  of  equation  (69)  are  the 
heat  energy  balance  at  the  soil  surface  and  specification  of  the  soil  temperature  at 
a  certain  depth.  The  heat  energy  balance  at  the  surface  is  given  by: 

R3  +  RlT-Rl1=H,  +  LE  +  G,  (70) 


where 

Rs  =  the  incoming  direct  solar  radiation  absorbed  by  the  surface, 
RlI  =  the  incoming  long-wave  radiation. 


and 


RlT  =  the  outgoing  long-wave  radiation. 


In  this  model,  the  loss  of  direct  solar  radiation  due  to  the  absorption  by  cloud  and 
the  small  effects  of  the  earth’s  annual  orbital  changes  are  taken  into  consideration, 
but  the  contribution  of  diffuse  solar  radiation  is  not.  The  sensible  heat  flux  H  j, 
latent  heat  flux  LE,  and  ground  heat  flux  G  5  are  given  by: 


(71) 


LE=  -p  Lu.Q, 


dT 

G=K. — 
dz_ 


(72) 

(73) 
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where 


p  =  the  air  density, 
u,  =  the  friction  velocity, 

T*  =  the  temperature  scale, 

Q.  =  the  water  vapor  scale, 

and 

subscript  G  =  the  value  at  the  ground  surface. 

Substituting  equation  (68)  to  (70)  into  (67),  we  obtain: 

R,  +  €R,  i  -€01*0  =  -p.C,u.T.-p.Lu.Q.  -  K.(aTyazJ  |  „  (74) 

where 

R^l=eor^  +  (\-e)R^l  (75) 


is  used, 

e  =  the  emissivity  of  the  surface, 
o  =  the  Stefan-Boltzman  constant. 


Garratt  and  Hicks  [12]  obtained  a  relationship  between  the  surface  temperature  and 
air  temperature  at  Zi  in  the  surface  layer: 


6(z,)-0^  z,+Zi.,  Zf. 

''  °  =-r[to-^  ■►In— 

r  k  Zn 


(76) 
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A  constant  value  of  0.1  m  is  assumed  for  Zq,  and  Zq,  is  obtained  from  a  relationship 
ln(zo/Zot)=2.  [12]  Using  equation  (56)  we  can  eliminate  T*  from  equation  (74)  to 
obtain: 


dT_ 


-€oV+m[0(z,)-r<;(/>y/>„r-]+A:/^)|^=O  (77) 


where 

h 

Pq  =  a  reference  pressure  (1000  mb) 

Pg  =  the  pressure  at  the  surface. 

Equation  (77)  can  be  linearized  by  noting  that 

Y  «+l  _  j*  n 

- ^«i 


(78) 


(79) 


where 

n  and  n+1  =  the  n  and  (n+l)***  time  steps  of  integration  (a  typical  time 
increment  used  in  integration  is  1  to  10  min). 

After  substitution  into  equation  (77)  of  the  approximation 

(V'/=4(r/)’v'-3(r/)^  (80) 
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we  obtain: 


G  s 


<^)Tr\'i)*R,*€R^\  +36eo(r<.")‘'-^m0"(z,)  (81) 


where  the  derivative  dT/^Zjlo  is  replaced  by  a  forward  finite-difference 
approximation, 

(1)  -  Tg"^')/AZs;  T5(1)  =  the  soil  temperature  at  the  first  grid  level  from 

the  surface; 


and 

AZj  =  the  distance  between  the  surface  and  the  first  grid  level  in  the  soil  layer. 

By  this  method,  equation  (69)  is  solved  numerically  in  finite-difference  form  by 
the  Laasonen  method.  [13]  By  this  method,  equation  (69)  reduces  to  AT  3=B 

where 

A  =  a  tridiagonal  matrix 
and 

B  =  a  column  vector. 

The  solution  is  conveniently  obtained  by  using  the  relation: 
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(T,),  =  +  F, 


(82) 


where 


(TJ,  =  the  soil  temperature  at  the  /*  grid  level  from  the  surface. 


Expressions  for  and  F j  for  /  >  1  are  determined  from  the  finite-difference  form 
of  equation  (69),  and  equation  (81)  determines  E,  and  F/.  From  equations  (81)  and 
(82),  we  obtain: 


Az„ 


i|  W\3  /  ^^0 

480(7(5  )"+/«(—)  ^ 

Pr^  AZ 


(83) 


and 


[R^-^ea(T^)Um0'(z^] 

4eo(TQf^p  +(^) 

Pn  Az. 


(84) 


Numerical  integration  of  equation  (69)  by  use  of  equations  (82)  through  (84)  is 
rapid  since  no  iteration  is  required. 

4.4  Radiational  Energy  Fluxes 


The  incoming  direct  solar  radiation  flux  to  an  inclined  surface  is  obtained  from 
Kondratyev:  [14] 
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(85) 


=Rf^A  +BcosQl + Cs/«Q]( 


where 

aVr^=  1.0001 10  +  0.03422  lcos(0o)  +  O.OO128Osin(0o)  +  0.0007 19cos(20o) 
+  O.OOOO77sin(20o)  (86) 

and 

0o=27ijy365  (87) 

A  =cosasin<E>sin6  +sina[cosT^(tan^>sin^sin6  -sin6secO)]  (88) 


J5=cosacos$cos6  +sinacos$sinT^cos6 


(89) 


and 


C=sinacos6sinY„ 

Yl 


(90) 


In  the  above  expressions, 

Ro  =  the  near  surface,  direct  solar  radiation  flux; 

Q  =  the  solar  hour  angle  (positive  clockwise  from  apparent  noon); 
=  the  latitude; 

6  =  the  declination  of  the  sun; 

a  =  the  angular  inclination  of  the  surface  to  the  horizontal  plane; 


53 


and 

Tn  =  the  azimuth  of  the  projection  of  the  normal  to  the  surface  on  the 
horizontal  plane,  as  counted  from  the  plane  of  the  meridian  (azimuth  is 
considered  positive  when  counted  clockwise). 

Because  the  maximum  change  in  the  solar  declination  6  in  24  h  is  less  than  0.5 
6  is  assumed  to  be  constant  during  a  given  day.  Paltridge  and  Platt  [15]  provide  a 
formula  to  compute  6  in  rads  by: 

6  =  0.006918  -  0.399912coseo+  0.070257sineo 

-  0.006758cos2eo  +  0.000907sin2eo 

-  0.002697cos3eo  +  0.001480sin3eo.  (91) 


Equation  (91)  estimates  6  with  a  maximum  error  of  0.0006  rd.  Solar  hour  Q  can 
be  obtained  if  the  longitude,  clocktime,  and  the  equation  of  time  are  known.  The 
equation  of  tune  is  the  difference  between  the  local  apparent  time  and  a  fixed  mean 
solar  time,  which  is  derived  from  the  motion  of  a  celestial  equation  at  a  rate  equal 
to  the  average  movement  of  the  sun.  The  solar  angle  Q  is  given  in  rads  by: 


Q= 


^(^,-12) 

12 


(92) 


where 

tj  =  the  true  solar  time  (local  apparent  time)  in  hours. 
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The  true  solar  time  is  obtained  from: 


t  =/  +A?/  +/ 

s  cA.  long  eq 


(93) 


where 

tc  Atjong,  and  the  clocktime,  the  longitude  correction,  and  the  equation  of 
time. 

The  longitude  correction  accounts  for  the  difference  between  the  local  meridian 
and  a  standard  meridian,  and  is  positive  if  the  local  meridian  is  east  of  standard. 
The  equation  of  time  is  as  follows:  [15] 

t  =il(0.000075 +0.001 868cos0„) 

71 

-O.O32O77sin0o-O.O14615co52eo-O.4O8495w20o)  ( 94) 


where 
U  =  hours 

00  =  defined  by  equation  (87). 

Equation  (94)  has  a  maximum  error,  compared  with  values  tabulated  in  the 
National  Almanac,  of  35  s  in  time.  The  amount  of  solar  radiation  reaching  the 
surface  is  much  less  than  that  at  the  top  of  the  atmosphere  due  to  many  factors, 
including  molecular  scattering  and  absorption  by  permanent  gases  such  as  oxygen, 
ozone,  and  carbon  dioxide.  The  effect  is  parameterized  by  Atwater  and  Brown 
[16],  who  modified  the  original  form  by  Kondratyev  to  include  the  effect  of  the 
forward  Rayleigh  scattering.  [17] 
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The  expression  is: 


(T=n4«s+0  ^i^[i  OAi  _n  i^(P-00Q949P+0.051  y,2^ 

cospo 


(95) 


where 

P  =  the  surface  pressure  in  megabytes. 

Other  important  factors  that  also  modify  the  amount  of  incoming  solar  radiation 
include  water  vapor,  clouds,  and  airborne  particles.  Parameterizations  for  clouds 
are  described  in  section  4.8.  Currently  Rq  is  calculated  from: 


(96) 


where 

Rfls  =  the  incoming  radiation  flux  at  the  top  of  the  atmosphere 

and  G  is  given  by  equation  (95).  The  zenith  angle  p  q  in  equation  (95)  is  determined 
from: 


cosp(j=sin$sin6+cos$cos6cosQ 

Finally,  Rj,  the  long-wave  incoming  radiation  at  the  surface,  is  computed 
according  to  the  following  formula: 

Rl-L  =  Roicosa  (98) 

where 

Ro^  =  the  long-wave  incoming  radiation  normal  to  horizontal  surface, 
a  =  the  angle  of  inclination  of  a  sloped  surface  given  by: 
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When  cloud  is  present  in  the  model  atmosphere,  the  surface  fluxes  of  both  short 
and  long  wave  radiations  are  modified  according  to  the  methods  described  in 
section  4.8. 

4.5  Upper  Boundary  Conditions 

Bovmdary  values  for  P,  U,  V,  0v,and  Q^,  along  the  upper  computational  boundary, 
are  given  by  the  values  yielded  from  the  analysis  of  the  GSM  plus  upper  air 
sounding  data.  The  values  are  updated  hourly  by  linearly  interpolating  between  the 
analysis  fields  at  Tq,  Tq  +  12  and  Tq  +  24.  Turbulence  and  vertical  velocity  are 
assumed  to  vanish  along  the  upper  boundary. 

4.6  Lateral  Boundary  Conditions 

The  lateral  boundary  values  for  U,  V,  0^,  Q^,  q^,  and  q  H  are  obtained  by 
integrating  the  corresponding  governing  equations  (35),  (36),  (42),  (43),  and  (45), 
except  that  variations  in  the  horizontal  directions  are  neglected.  Variables  U,  V, 
0v>  Qw»  smoothed  at  each  time  step  by  using  the  values  at  four 

neighboring  points,  such  as.: 

<Pv=(l-^)‘P//0.25A((p.,j>(p..,^.+(p.^^  (100) 

where 

(p  =  either  U,  V,  0^,  q^  q^/, 

A  =  smoothing  parameter  (=0.25). 


A  similar  expression,  using  only  three  neighboring  points,  is  applied  to  the  values 
at  the  lateral  boundaries.  In  addition,  changes  have  been  made  in  the  code  that 
damp  the  effects  of  terrain  slopes  at  the  boundaries.  This  makes  it  possible  to  put 
the  boundaries  in  very  complex  terrain  without  constructing  an  artificial  apron 
around  the  model  domain.  Also,  the  model  code  forces  the  first  grid  vertical  cell 
to  be  at  4  m  AGL.  This  is  to  ensure  that  the  lowest  layer  is  sufficiently  near  the 
ground  to  generate  appropriate  slope  flows. 

4.7  Diagnostic  Pressure  Distribution 

In  the  original  code  of  HOTMAC,  three-dimensional  pressure  distribution  is 
diagnostically  calculated  at  each  time  step  from  the  initial  distribution  of  pressure 
at  the  top  of  the  model  atmosphere  and  the  three-dimensional  temperature 
distribution  that  is  calculated  prognostically.  However,  we  found  that  by  this 
method  the  pressure  distribution  sometimes  becomes  unrealistic. 

By  using  a  time-dependent  distribution  of  pressure  at  the  top  of  the  model 
atmosphere,  as  described  in  the  following,  a  more  dependable  pressure  distribution 
is  obtained. 


Initially,  the  surface  pressure  distribution  is  calculated  from  the  GSM  and  upper 
air  sounding  data.  Using  this  data,  the  three-dimensional  pressure  distribution  over 
the  model  domain  is  calculated  as  follows.  A  modified  pressure  is  defined  as: 


n=ci— f 

^0 


(101) 
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where 


Pq  =  a  reference  pressure  (1000  mb), 

K  =  the  ratio, 

Cp  =  the  specific  heat  of  the  dry  air  at  constant  pressure, 

Rj  =  gas  constant  of  dry  air. 

The  hydrostatic  equation  in  Cartesian  coordinates  can  be  written  as: 

5n  _  g 

dz  0^ 

In  the  transformed  coordinate, 

an  _an  dz  _  (h-z^) 

dz*  dz  dz*  ^  HQ^ 


(102) 


(103) 


From  equation  (103),  the  modified  pressure  at  the  top  of  model  atmosphere  is 
calculated  as: 


n(^=n(0)- 


g(H-Zg)  J^dz 
H  *'0  0v 


(104) 


where  the  modified  pressure  at  the  surface  n(0)  is  calculated  from  a  surface 
pressure  value  as: 

(105) 
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From  the  initial  and  boundary  condition  analyzed  fields  of  the  surface  pressure 
from  the  GSM  (and  at  initial  time,  the  upper  air  sounding  data  as  well),  the 
modified  pressure  distribution  at  the  top  of  the  model  atmosphere  is  calculated  by 
using  equation  (101)  at  the  initial ,  12-,  and  24-h  forecast  time.  The  modified 
pressure  at  the  top  of  the  model  atmosphere  at  time  t  after  the  initialization  for  t 
less  than  12  h  is  now  calculated  as: 


n(/0rn(i^)o 


n(/o,2-n(/^o 

+ - if - 

12 


(106) 


where  the  subscripts  0  and  12  stand  for  t=0  and  12  hours,  respectively.  A  similar 
equation  is  used  for  t  between  12  and  24  hours. 


The  modified  pressure  at  the  level  z*  at  time  t  is  calculated  as: 

n(z  *)=n(m+-^^^ 

H  0v 


(107) 


4.8  Radiative  Effects  in  Clouds 

Solar  and  longwave  radiation  play  important  roles  in  the  formation  and  dissipation 
of  clouds  and  fog.  For  example,  radiation  fog  forms  when  moist  air  near  the 
surface  condenses  as  the  ground  temperature  decreases  due  to  longwave 
radiational  cooling.  As  the  sun  rises,  fog  reflects,  absorbs,  and  transfers  shortwave 
radiation.  The  solar  energy  absorbed  in  a  fog  layer  heats  the  air  and  converts  liquid 
water  to  water  vapor.  The  solar  energy  transmitted  through  fog  heats  the  ground 
and  increases  ground  temperature.  The  warmed  ground,  in  turn,  heats  the  air  above 
and  dissipates  the  fog.  Interaction  between  radiative  transfer  and  cloud 
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development  is  understood  qualitatively.  However,  a  quantitative  description  is 
complex  and  involves  considerable  computations.  [18]  [19] 

Hanson  and  Derr  proposed  a  parameterized  solar  radiation  scheme  whose 
parameters  were  obtained  by  curve-fitting  to  numerical  radiative-transfer  results 
using  the  atmospheric  radiation  (ATRAD)  narrow  band  model.  [20],  [21]  This 
parameterized  method  is  simple,  yet  reproduced  solar  flux  profiles  within  a  single 
cloud  layer  that  were  in  good  agreement  with  the  numerical  results  obtained  from 
ATRAD. 


Motivated  by  the  success  of  the  solar  radiation  scheme,  Hanson  and  Derr  proposed 
an  infrared  radiation  scheme  expressed  by  exponential  functions  whose  decay 
parameters  were  determined  by  the  emissive  method.  [22]  This  method  reproduced 
infrared  flux  profiles  within  layered  clouds  that  were  in  good  agreement  with  the 
numerical  results  and  observations  for  thick  clouds  ( about  800  m).  However,  the 
parameterization  overestimated  the  flux  decrease  at  the  cloud  top  and 
underestimated  the  cloud  base  warming  compared  to  the  numerical  model  for  thin 
clouds  (about  300  m). 

The  BFM  has  adopted  the  Hanson  and  Derr  parameterization  scheme  because  of 
simplicity  and  because  it  can  produce  flux  profiles  that  are  in  good  agreement  with 
observations  and  numerical  results.  [23] 

4.8.1  Shortwave  Radiation 


Following  the  two-stream  model  solutions  by  Stephens  et  al.,  reflection  (R,.), 
transmission  (T^)  and  absorption  (A^)  are  expressed  as:  [24] 

(1)  Nonabsorbing  medium  [  t3o=  1  (A  ^  0.75  pm )] 


RM 


1  +P(^o)VlAo 


(108) 
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TXHo)  =  1  -  Re(^o)  (109) 

(2)  Absorbing  medium  [  GJo  ^  1  (A  ^  0.75  ^im  )] 

Re(^o)  =  [(u^  - 1  )exp(T,ff)  -  exp(-T,ff)]/R  (110) 

T,(Po)  =  4u/R  (111) 

Ab(Po)  =  1  -  Re(lio)  -  Tr(Po)  (1 12) 

where 

=  [1  -  ©0+  2p(iio)ao]/(l  -  ©o)  (1 13) 

{(1  -  ©o)[l  -  ao  +  2P(iio)t3o]}'%/po  (114) 

R  =  (u  +  l)"exp(T,ff)  -  (u  - 1  )"exp(-T,ff)  .  (115) 

In  the  above  expression, 

Po  =  the  zenith  angle, 

=  the  optical  thickness  of  the  cloud, 
tUg  =  the  single-scattering  albedo. 


P  =  the  backscattered  fraction  of  monodirectional  incident  radiation. 

The  values  for  ©g  and  P  are  tabulated  as  functions  of  Tg  and  pg. 

The  solar  radiation  flux  profile  within  the  layered  cloud  is  computed  from  the 
following  equations:  [20] 
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F(z)=Fb-  (Fb-  Fc)[1  -  exp(-(ZB  -  z)/As)]/ 


{1  -  exp[  -  (Zb  -  Zc)/XJ}  (1 16) 

where 

As  =  a(|io)W  +  b(^o)(  1  -  exp{  -[yW  +  c(po)])  (117) 

Fc^FB+AbFei  (118) 

a(po)  =  -0.022  +  0.038(1  -  Po)  (1 19) 

b(po)=  56.8- 14.7(1 -po)  (120) 

c(po)=  1.07 -1.15(1 -po)  (121) 


In  the  above  expressions  W  is  the  total  cloud  liquid  water  content,  which  is  given 
as: 


p(z)Q/z)dz 


(122) 


where 


p(z)  =  the  air  density, 

Zg  =  the  cloud  top, 
z^  =  the  cloud  base. 

The  units  of  W  are  gm‘^  and  y=0.021. 

Fb=  -F3i{l.R(po)} 


(123) 
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where 


Fgi  =  the  downward  solar  flux  at  the  top  of  the  cloud. 

The  derivation  of  cloud  water  mixing  ratio  Q^,  is  described  in  the  appendix. 


4.8.2  Longwave  Radiation 

The  parameterization  of  longwave  radiation  relies  on  the  values  of  external 
conditions:  Tg  and  Tc  are  the  absolute  temperatures  at  the  cloud  top  and  the  cloud 
base,  respectively.  Ggi  is  the  downward  longwave  flux  at  the  cloud  top. 
GqT  is  the  upward  longwave  flux  at  the  cloud  base.  The  infrared  flux  profiles  in  the 
layered  cloud  are  given  as: 

G(z)  =  GLexpf  -  (z  -  zJ/A.  J  +  Gycxpf  -  (zg  -  z)/Au]  (124) 


where 


Gl  =  {  Go  -  G,exp[  -  (Zb-  z,)/Au]}/D 

(125) 

Gu  =  {  G,  -  Goexpf  -  (Zg  -  z,)/AJ}/D 

(126) 

D  =  1  -exp{-(zB-zJ/AN} 

(127) 

J-=-L+  ^ 

A«,  A.,  A, 

(128) 
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^0  (®c  ■ 

(129) 

Gi  =  (GJ  +  Bb  -  2Bc)exp(  -  T|T)  +  Bg  -  Ggi 

(130) 

TiT,i  =  aT,iW 

(131) 

at  =  0.13  and  ai  =  0.158 

(132) 

Be  =  oTc^ 

(133) 

Bb  ~  oTb 

(134) 

Xl  =  70W/(W-W'^  +  2.67) 

(135) 

Xu  =  140W-°'^ 

(136) 

where 

o  =  the  Stephan-Bolt2aiiann  constant. 

4.9  Cloud  Condensation 

Conservation  equations  for  the  mixing  ratios  of  water  vapor  and  cloud  water,  and  the 
potential  temperature,  are  written  as  follows; 

dO 

-^=-(Cond)^  (137) 

dt 


65 


(138) 


dt 


={Cond)^ 


00 

dt 


—iCond) 

^  s 


V 


(139) 


where 

Qy  Qc  =  the  mixing  ratios  of  water  vapor  and  cloud  water. 

0  =  potential  temperature, 

T  =  absolute  temperature, 

Ly  =  latent  heat  of  condensation, 

Cp  =  specific  heat  of  dry  air  at  constant  pressure, 

and 

(Cond)y  =  the  rate  change  of  the  mixing  ratio  due  to  condensation. 


Define: 


Qw-Qv  +  Qc 


L 

0=e-^^Q 

T  c 


(140) 

(141) ) 
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where 


=  the  total  water 

0,  =  liquid  water  potential  temperature.  [25] 

In  model  operation,  is  initially  zero,  and  comes  from  initial  input  data.  By  adding 
equations  (137)  and  (138),  we  obtain: 

dO 

— =0  (142) 

dt 

By  subtracting  equation  (138)  x  (6/T)(LyCp)Q^from  equation  (139),  we  obtain: 

30, 

— ^=0  (143) 

dt 


Thus,  and  0,  conserve  their  quantity  even  when  a  phase  change  occurs. 

In  order  to  recover  the  potential  (or  absolute)  temperature  and  the  mixing  ratios  of 
water  vapor  (QJ  and  cloud  liquid  water  (Q^),  the  probabilty  density  function,  G, 
proposed  by  Sommeria  and  Deardoff  [26]  and  Mellor,  [27]  is  used.  The  probability 
density  function  in  assumed  to  be  Gaussian  such  that: 


G=- 


1 


2\\I2 


exp[- 


1 


(l-r2)  2ae/ 


)] 


(144) 


67 


where 


0/  -  fluctuation  of  liquid  water  potential  temperature, 
the  fluctuation  of  total  water  mixing  ratio, 

and 

(145) 


(146) 


Go, a 

0/  qw 


(147) 


The  local  condensation  is  assumed  to  be  given  by: 

Q/  =  (Qw  -  Qs)H(Qw  -  Qs)  (148) 

where 

Qs  =  saturation  mixing  ratio, 

H(x)  =  Heaviside  function. 
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defined  as: 


H(x)=0,  x<0 
H(x)=  1,  x>0 

The  following  parameters  are  defined  as: 


(149) 


(150) 


b=a 


<T> 

<0> 


sl,T 


(151) 


AQ  =  Qw-Qsi 


L  Q, 

Q,,t=(—)t  r  =0.622— ^ 

^sl,T  y  gj,JT-T,  ^  2 

d  ■*  / 


Q,=  0.622e3(T,)/(P-e3(T,)) 


e^(7)=6.11exp[ 


—(—-—)] 
273  T, 


and 

T,  -  (P/P/0, 


(152) 

(153) 

(154) 

(155) 

(156) 
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where 


Qsi  =  mean  saturation  mixing  ratio  of  water  vapor  at  T, ; 

0/  =  liquid  water  potential  temperature  defined  by  equation  (138); 

es(T,)  =  saturation  water  vapor  pressure  at  T/  ; 


and 

=  the  gas  constant  for  water  vapor. 

A  function  R  which  indicates  a  fraction  of  cloud  coverage  for  a  given  volume  of  air  is 
given  by  Mellor;  [27] 


Jl=f  7  ~H(,Q.-Q,)GdQ^deri[UerAQfl^)] 

J  —oo  J  — oo  ? 

erf{x)  =-^  f  "'expC  -y  ^)dy 


(157) 

(158) 


Cloud  water  mixing  ratio  is  given  by: 

+-^exp(  -^)] 
2 


(159) 


where 


Q\=a 


A0 


2a 


S 


(160) 
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(161) 


o/  =^{a  V  -^abq^e,+b  ^0 /) 


Furthermore, 

(162) 


and 

ufljaujq^ -bUjQj  =R'  (1 63) 


where 

J!'-/e-^-^exp[-^].  (164) 

20,^2^  2 


However,  due  apparently  in  part  to  coarse  vertical  grid  spacing,  and  the  nudging 
discussed  later  in  section  5.2  the  values  of  o^  calculated  by  equation  (161)  are  typically 
too  low.  In  the  current  BFM,  a  constant  value  of  0.000275  is  used  instead.  This  was 
found  to  be  typical  of  partly  cloudy  skies  by  Sommeria  and  Deardoff.  [26]  If  model 
nudging  is  turned  off,  equation  (161)  seems  to  yield  more  realistic  values.  Finally,  q 
the  variance  of  the  cloud  water  mixing  ratio,  is  expressed  as: 

ql  Q  Q  \  Qf 

(165) 

4o^  2a^  2o^  ^  2 
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Figure  7  shows  R,  Q^20j,  and  as  a  function  of  Q|  obtained  according  to 

equations  (157),  (159),  and  (165),  respectively.  As  seen  from  Figure  7,  a  fraction  of 
cloud  coverage  R  varies  gradually  with  Q]  and  takes  non-zero  values  even  if  Q]  is 
negative.  In  other  words,  clouds  can  exist  even  if  the  mixing  ratio  of  water  vapor 
averaged  over  a  grid  volume  is  not  saturated.  This  is  realistic  because  the  grid  spacing 
normally  used  in  mesoscale  models  is  larger  than  the  size  of  small  clouds.  Therefore, 
the  present  cloud  model  can  use  a  relatively  large  grid  spacing  that  could  save 
computational  time  substantially.  A  statistical  cloud  model  such  as  the  present  one 
avoids  the  ambiguous  condensation  criteria  often  used  by  coarse  grid  models  in  which 
saturation  values  are  lowered  arbitrarily  to  compensate  for  the  amount  of  cloud  that  is 
not  resolved  by  the  grid. 


Q|  =  0AQ/(2crs) 


Figure  7.  Relationships  showing  R,  Q  ^I2a„  and  as  a  function  of 

Q.  (=  a  (Qw-Qsi)  /2o,). 
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In  the  BFM,  the  total  cloud  water  content  (W)  from  the  ground  surface  to  the  top  of  the 
model  atmosphere  is  calculated  and  the  horizontal  distribution  of  W  is  displayed.  W 
is  calculated  as: 


H+z  -z 

gmax 


Q 

u  Jo 


.dz 


H 


(166) 


where 

p  =  air  density, 

Zgmax  =  highest  terrain  elevation  value  in  the  BFM  5 1  by  5 1  model  domain. 

4.10  Non-Convective  Precipitation  Rate 

Because  microphysical  processes  of  non-convective  precipitation  formation  are  not 
included  in  the  model,  the  non-convective  precipitation  rate  is  parameterized  as  a 
function  of  cloud  liquid  water.  The  scheme  developed  by  Sundqvist  et  al.  for  non- 
convective  precipitation  rate  is  incorporated  into  the  BFM.  [28]  The  basic  assumption 
is  that  the  denser  the  cloud,  the  higher  the  precipitation  rate. 


The  rate  of  release  of  non-convective  precipitation  is  described  by 


P=C„0[l-exp(-(-^)'] 


(167) 


c,cr 
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where 


I/Cq  =  characteristic  time  for  the  conversion  of  cloud  droplets  into  precipitation 
particles  (raindrops,  or  snow  particles), 

R  =  fraction  of  cloud  coverage  described  in  the  previous  section, 

Qc  „  =  threshold  value  for  cloud  water,  which  Q  ^/R  must  exceed  before  the 
release  of  precipitation  can  become  efficient. 

The  parameter  should  have  a  value  typical  of  individual  cloud  type,  which  is 
invariant  to  grid  resolution. 

The  rate  of  non-convective  precipitation,  at  a  z*-level  is  given  by 
^  u  H-^z  -z  T- 

"iPdz  ■  (168) 

where 

H  =  depth  of  the  model  atmosphere, 

Zgmax  =  highest  terrain  elevation  in  the  BFM  model  5 1  by  5 1  domain, 

Zg  =  terrain  elevation,  and  p  is  the  air  density. 

To  simulate  the  coalescence  process,  Sundqvist  et  al.  introduced  an  additional 
parameter  F„.  [28]  This  parameter  increases  with  the  rate  of  precipitation  and 
multiplies  Cq  and  divides  Q^„.  The  relation  is: 

F„(z*)=l+C,Pr(z7^  •  (169) 
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Similarly,  Cq  increases  and  decreases  by  a  temperature  function  Fbf  when  the 
temperature  is  lower  than  -5  °C  in  clouds  containing  a  mixture  of  droplets  and  ice 
crystals  (Bergeron-Findeisen  mechanism).  This  formulation  is: 

Fbf=  1  +  €2(268.0  -  T(z*))*^  .  (170) 

The  two  modified  parameters  Cqf  and  Q(c,cr)F  become,  respectively: 

Cof  =  CooF„=Fbf  (171) 

and 

Q(c.cr)F“Qc,c/(F~'FBF)  •  (172) 

The  effects  of  evaporation  as  precipitation  falls  through  subsaturated  layers  is  not 
included.  Also,  the  non-convective  precipitation  rate  at  the  surface  is  calculated  as: 

p  (0)=(:^!^£l:!£)±  *)CoFe,(^  *)[1  -exp(-(-^-^)''')]d!^  •  (173) 

H  ^Qic,cr)F 


where 

=  density  of  water. 

Pr(0)  is  calculated  in  the  unit  of  (mm/h)  or  (in/h). 
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4.11  Sea  Surface  Temperature  and  Soil  Characteristics 


4.11.1  Sea  Surface  Temperature 

In  the  BFM,  water  temperature  stays  independent  of  diurnal  change  of  the  solar 
angle,  but  the  monthly  average  global  sea  surface  temperature  distribution  data, 
obtained  from  the  U.S.  Air  Force  Combat  Climatology  Center,  is  used  to  give  the 
temperatures  of  sea  water.  The  data  is  given  at  about  every  0.5°.  The  temperatures 
of  the  grids  over  sea  surface  in  the  model  domain  are  given  by  those  adjacent  in 
the  data  set.  Currently,  water  temperatures  of  inland  water  bodies  such  as  the  Great 
Lakes  and  the  Salt  Lake  must  be  given  by  changing  the  source  program  before 
application  of  the  model. 

4.11.2  Surface  Characteristics 

Ground  surface  albedo  and  emissivity,  specific  heat,  density,  and  heat  conductivity 
of  soil  are  given  as  a  fimction  of  the  month  of  year,  and  latitudes  of  the  model 
grids.  The  values  of  these  parameters  are  roughly  estimated  from  the  book  by 
Pielke.  [28] 

For  snow  cases,  the  user  can  specify  that  snow  can  be  assumed  over  the  entire 
grids  above  a  particular  altitude  by  manually  setting  flags  in  the  file  snow.d.  For 
the  land  grid  points  that  are  snow  covered,  albedos,  emissivities  and  soil  heat 
conductivities  are  internally  hardcoded  to  those  typical  of  snow,  and  soil 
temperature  is  fixed  at  -  1  °C.  This  has  not  yet  been  added  as  part  of  the  IMETS 
BFM  execution  GUI,  so  for  now  this  file  is  set  to  “no  snow.” 
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5.  Initialization  and  Assimilation  of  Data 


5.1  Initialization 

Initialization  of  the  BFM  is  done  as  follows.  The  initial  wind  field  over  the  entire 
model  domain  is  given  uniformly  by  the  average  wind  vector  field  calculated  by 
the  mean  wind  direction  at  the  tenth  layer  (Z*=  1,195  m),  and  the  maximum  wind 
speed  of  the  following  two: 

Mean  wind  speed  at  the  tenth  layer,  and 
(0.4/k)log([z‘  +  0.1}/0.1) 

where 

k  =  Karman  constant. 

The  average  wind  direction  and  speed  at  the  tenth  layer  are  calculated  fi’om  the 
three-dimensional  data  obtained  fi’om  the  analysis  of  the  GSM,  and  upper  air 
sounding  data.  Initial  fields  of  potential  temperature  and  water  vapor  mixing  ratio 
are  also  obtained  fi'om  the  same  data. 

Suppose  a  forecast  calculation  is  initialized  at  time  Tq.  Pre-calculation  will  start  at 
time  To-3  h,  and  for  3  h  from  Tq  -  3  to  Tothe  model  fields  are  dynamically  adjusted 
to  the  initial  fields  by  the  nudging  method  which  is  described  in  the  following 
sections.  During  the  3-h  initialization  period,  surface  data  observed  at  the  time  Tq 
is  also  nudged  to  the  model  fields.  The  nudging  method  of  surface  data  will  be 
described  later. 

When  the  GSM  forecast  data  for  time  Tq,  Tq  +  12,  and  Tq  +  24  h  are  available,  the 
hourly  lateral  boundary  condition  data  between  T  q  and  Tq  +  12,  and/or  Tq  +  12  and 
Tq  +  24  are  calculated  by  a  linear  interpolation  method.  At  time  Tq,  the  data  for  Tq 
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+  1  are  assimilated  in  for  1  h,  and  this  process  is  repeated  for  an  entire  forecast 
period. 

When  only  upper  air  sounding  and  surface  data  are  available,  initialization  is  done 
in  a  similar  manner  as  when  the  GSM  data  is  available,  but  after  the  time  Tq,  the 
magnitudes  of  nudging  parameters  C„  are  reduced  gradually  as  exp(-a't),  where  t 
is  the  forecast  duration  after  Tq,  and  a  is  a  coefficient  determined  empirically.  The 
forecast  period  is  limited  to  6  h. 

5.2  Nudging  Method 

To  assimilate  temporal  and  spatial  changes  of  synoptic  flows,  which  are  obtained 
from  the  GSM  and  upper  air  sounding  data,  nudging  terms  are  added  to  the  model 
equations  (35),  (36),  (45),  and  (46)  at  all  grid  points,  although  the  model  can  be 
run  with  nudging  only  at  the  lateral  boundaries  as  well.  By  adding  the  nudging 
terms,  the  equations  of  horizontal  wind  vector  components  are  modified  as: 

^‘Fx*CJ.V,-U)  (174) 


(175) 

where 

C„  =  nudging  coefficient, 

U,  and  V,  =  target  wind  components  of  (respectively)  x-  and  y-  directions. 
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Target  winds  are  described  in  the  next  section.  F ,  and  Fj  are  the  right-hand  side 
terms  of  equations  (35)  and  (36),  respectively.  The  equations  of  potential 
temperature  deviation  and  mixing  ratio  are  written  with  nudging  terms  as  follows: 

-^=F3.C„(60,,,,-60,)  (176) 


dt 


^A'^^n^Qv,obs  Q-J 


(177) 


where 

60y  =  deviation  of  virtual  potential  temperature  from  the  base  state  virtual 

potential  temperature, 

<0y>  =  defined  earlier, 

Qy  =  water  vapor  mixing  ratio, 

obs  =  observed  values,  which  are  obtained  from  three-dimensional 

analysis  of  the  GSM  and  upper  air  soimding  data, 

F3  =  right-hand  side  terms  of  equation  (45), 

F4  =  right-hand  side  of  (46). 

Recall  that  in  the  BFM,  <0^>  is  the  initially  analyzed  three-dimensional  virtual 
potential  temperature  field,  rather  than  the  initial  horizontally  averaged  values 
used  by  modelers  such  as  Yamada.  [5] 

In  the  current  version  of  the  BFM,  the  nudging  are  enforced  to  the  vertical  layers 
above  z*  =  150  m  for  potential  temperature  and  moisture,  and  above  z  *  14  m  for 

wind  components.  These  are  determined  based  on  our  long  term  experiences  with 
nudging.  The  value  of  the  nudging  coefficient  is  set  to  0.0003/s,  roughly  the 
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magnitude  of  the  coriolis  force.  This  is  a  value  that  has  typically  been  used  by 
many  modelers  such  as  Pielke,  Yamada,  and  Zack  et  al.  [29],  [4],  [30]  This  value 
seems  to  properly  assimilate  the  large  scale  perturbations  from  the  GSM,  while  at 
the  same  time  allows  the  surface  physics  of  the  mesoscale  model  to  have 
significant  influence  within  the  boundary  layer,  particularly  near  the  surface. 

5.3  Target  Winds 

Comparisons  in  the  simulated  wind  flelds,  found  by  nudging  to  target  wind 
components  Uj  and  ,and  by  nudging  to  the  observed  wind  components  U  ^  and 
^obs  were  conducted.  [31]  It  showed  that  nudging  to  the  target  wind  components 
produces  better  agreement  between  simulated  and  observed  upper  winds  than 
nudging  to  the  observed  wind  components. 

Target  wind  components  U,  and  Vj  are  derived  as  follows.  The  equations  of 
motion  for  horizontal  wind  vector  components  with  the  target  winds  under  no 
frictional  force  can  be  written  as: 

^=AV-V^hC„(U-U)  (178) 


^=-AU-u^)^c„{v-v) 


(179) 


The  solutions  to  the  above  equations  at  equilibrium  (t  -*<»)  are  given  by  the 
following: 


t/= 


(180) 
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(181) 


V= 


cl-f  . 


Letting  U  and  V  approach  U^bs  and  V^bs  (respectively)  and  solving  for  U,  and 
the  following  equations  are  obtained: 


(182) 

(183) 

where 

Uobs  =  the  observed  wind  components  in  the  x-  direction 
Vobs  =  the  observed  wind  components  in  the  y-  direction 

and 

Ug  =  geostrophic  wind  components  in  the  x-  direction 
Vg  =  geostrophic  wind  components  in  the  y-  directions. 

Uj  and  Vj  are  generally  different  from  the  corresponding  large-scale  wind 
components.  Observed  winds  might  be  used  as  target  winds  if  the  Coriolis  force 
were  absent,  or  if  the  observed  winds  are  identical  to  the  geostrophic  winds.  If 
only  the  observed  winds  are  used  in  the  nudging  in  all  other  cases,  the  solutions 
will  generally  be  different  from  the  observations. 
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The  physical  meaning  of  the  target  winds  is  that  the  solutions  of  the  equations  of 
motion,  with  the  target  winds,  become  identical  to  the  observed  winds  in  the 
absence  of  frictional  effects.  Thus,  modeled  winds  should  converge  to 
observations  in  the  layers  far  above  the  boundaiy  layer  where  frictional  effects  are 
iicgligible.  On  the  other  hand,  atmospheric  turbulence  in  the  boundary  layer  is 
significant  because  of  frictional  effects.  The  nudging  terms  play  relatively  minor 
roles.  In  summary,  the  nudging  terms  enforce  the  model  winds  to  match 
observations  in  the  free  atmosphere,  but  they  play  a  relatively  minor  role  in  the 
boundary  layer. 


5.4  Surface  Data  Nudging 


Individual  surface  data  within  the  model  domain  observed  at  initialization  time  Tq 
is  assimilated  into  model  calculation  at  the  grid  points  adjacent  to  the  surface 
observation. 


The  &ird  and  forth  layer  data  (z*—  6  and  10  m,  respectively)  of  the 
three-dimensional  data  for  initialization  calculated  from  the  GSM  and/or  upper 
air  sounding  data  at  time  Tq  are  modified  as  follows: 


^  ill  ^ 


1 


/=!  rf  M  rf 


(184) 


where 

'l^new  ~  I16W  values  of  a  parameter  iji, 
iji/  =  iji  at  a  surface  station  /. 
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Vj  =  distance  from  /  to  a  grid  point  (i,  j). 


In  order  to  limit  spatially  the  influences  of  surface  data,  the  nudging  parameter  at 
a  grid  point  (i,  j)  is  calculated  as: 


1=1 


forr/<R  (185) 


and 


C 


„.„ew(ij)=0.0 


forr/^i?  (186) 


where 

C,  =  empirical  parameter  (0.1). 

The  critical  distance  R  is  40  km  for  wind  and  20  km  for  potential  temperature 
deviation  and  water  vapor  mixing  ratio.  These  values  are  determined  empirically. 
The  nudging  parameter  C ,  is  about  three  orders  of  magnitude  greater  than  C  „  to 
emphasize  the  effects  on  the  model  fields  of  surface  data. 

Surface  data  nudging  is  done  for  3  h  from  Tq  -  3  to  T  q.  After  T  q,  instead  of  a 
sudden  assignment  of  zero  values  to  the  influences  of  surface  data  are 
gradually  decreased  by  multiplying  exp(-kt)  to  C^.  Here  k  is  an  empirical 
coefficient,  and  t  is  the  time  after  Tq. 
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6.  Initialization  Data  and  Some  Examples  of  the  BFM 


In  this  section,  the  qualities  of  the  GSM  forecast  data  used  as  the  BFM  lateral 
boundaiy  conditions,  the  initialization  data  produced  by  the  three-dimensional 
analysis,  and  some  examples  of  the  BFM  output  are  shown.  Preliminaiy  studies 
that  attempt  to  make  a  statistical  evaluation  of  the  BFM  output  have  been  done 
by  Henmi  et  al.  [32]  and  Knapp  and  Dumais.  [33]  Henmi  et  al.  studied  the 
performance  of  the  BFM  using  meteorological  data  observed  at  White  Sands 
Missile  Range  (WSMR),  NM,  and  showed  that  incorporation  of  surface  data  into 
the  initial  fields  seems  to  produce  more  accurate  BFM  forecast  results.  Knapp  and 
Dumais  made  comparison  studies  of  observed  data  including  upper  air  and  surface 
data  with  the  BFM,  and  the  GSM  output  over  a  250  x  250  km,  and  a  500  x  500 
km  area  centered  at  Colorado  Springs,  CO.  Statistical  parameters  such  as  the 
correlation  coefficients  and  mean  absolute  residuals  between  observed  data  and 
model  output  data  are  calculated  for  horizontal  wind  vector  components, 
temperature,  and  dew-point  temperature.  Knapp  and  Dumais  concluded  that  the 
BFM  forecasts  of  wind,  temperature  and  dew  points  at  the  surface  (and  to  a  lesser 
degree  with  height  above  the  surface)  are  significantly  better  than  those  by  the 
GSM  forecasts  being  used  as  the  BFM  lateral  boundary  conditions,  particularly 
over  complex  terrain.  As  to  the  effects  of  grid  spacing  of  the  BFM,  they  did  not 
find  any  significant  improvement  by  adapting  the  2.5  or  5  km  grid  spacing  in  place 
of  the  10  km.  It  might  be  expected  that  these  results  could  change  when  a 
convective  parameterization  scheme  is  incorporated.  This  problem  will  be  further 
studied  in  the  future.  Currently,  extensive  studies  using  the  data  archived  over 
different  geographical  and  climatological  regions  (Bosnia,  Korea,  Florida, 
southern  California,  and  the  northeast  U.S.)  during  different  seasons  of  the  year 
are  under  way  and  the  results  will  be  presented  in  the  future. 
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6.1  Qualities  of  the  GSM  Forecast  Data 


Initialization  and  time-dependent  boundary  condition  data  for  the  BFM  are 
provided  by  the  GSM  forecast  data.  Because  the  GSM  analysis  and  forecast  data 
are  delivered  about  6  to  7  h  after  the  forecast  base  time,  the  GSM  analysis  data 
fields  are  not  used  for  operational  initialization  of  the  BFM  forecast  calculations; 
the  12,  24,  36  and  48-h  GSM  forecast  data  fields  are  used.  For  example,  to 
initialize  a  BFM  forecast  at  12  GMT,  the  12  GMT  upper  air  radiosonde  and 
surface  sensor  data,  along  with  the  GSM  gridded  12-h  forecast  fields  fi'om  the 
previous  00  GMT  initialized  run  are  used.  Time-dependent  boundary  values  for 
the  12  GMT  run  come  fi'om  the  GSM  00  GMT  base  time  24, 36  and  48-h  gridded 
forecast  fields. 

In  this  section,  the  12  and  24-h  forecast  data  generated  by  the  GSM  have  been 
spatially  interpolated  to  the  radiosonde  sounding  site  locations,  and  compared 
with  the  observation.  The  GSM  data  calculated  at  different  pressure  levels  (200, 
300,  500, 700,  and  850  mb)  are  compared  with  the  corresponding  observed  data. 
The  data  at  1000  mb  are  excluded  due  to  the  lack  of  sounding  data  at  this  level  for 
the  stations  in  the  areas  selected  for  study.  The  data  are  obtained  fi-om  the  1,600 
X  1,600  km  areas  centered  at  the  following  geographical  areas  and  dates/hours  as 
shown  in  table  2. 
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Table  2.  Locations,  dates,  and  hours  of  upper  air  sounding  data  used  for 
comparison  with  GSM  data _ 


Location 

Longitude 

Latitude 

Base  Time+ 

Dates/Hour 

Tuzla,  Bosina 

18.657  E 

44.470  N 

12/13/12Z,12/14/00 

Z12/14/12Z,12/15/0 

oz 

12/19/12Z,12/20/00 

Z 

12/20/12Z,12/21/00 

Z 

12/21/12Z 

Denver,  CO 

104.717W 

38.17N 

7/24/12Z,7/31/12Z 

8/l/00Z,8/l/12Z 

8/2/OOZ,  8/7/1 2Z 

8/8/OOZ 

Ft.  Irwin,  CA 

116.066W 

34.257  N 

l/4/12Z,l/5/00Z 

1/10/12Z,1/11/00Z 

l/23/00Z,l/24/12Z 

1/25/OOZ,  1/31/12Z 

In  table  2,  12/13/12Z  means  that  the  GSM  base  time  is  13  December,  12  GMT. 
The  data  is  taken  during  the  periods  of  July  1995  through  January  1996. 

The  correlation  coefficients  r,  and  linear  relationships  y  =  A  +  Bx  between  the 
GSM  data  (y)  and  the  upper  air  sounding  data  (x)  at  different  pressure  levels,  and 
those  of  the  entire  data  set  are  calculated  for  the  12  and  24-  hour  forecast  data. 


87 


The  results  are  given  in  table  3.  Scatter  diagrams  showing  the  relationships 
between  the  GSM  data  and  upper  air  sounding  data  for  the  entire  data  sets  are 
shown  in  figure  8A  through  8D. 
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Table  3.  Correlation  coefficients  between  the  GSM  data  and  upper  air  data 


Pressure  levels 

(mb) 

Temperature 

12  h  24  h 

Dew  pt.  Temp. 

12  h  24  h 

x_comp.  of 

wind,  u 

12h  24h 

y-comp.  of 

wind,  V 

12  h  24  h 

850 

r=.83  r=.85 

A=.71  A=.58 

B=.86  B=.84 

r=.66  r=.59 

A=.71  A=-.28 

B=71  B=63 

r=.48  r=.52 

A=.23  A=1.03 

B=.58  B=65 

r=.28  r=.38 

A=-.36  A=-1.97 

B=.28  B=.40 

700 

r=.78  r=.82 

A=-1.5  A=-1.2 

B=.86  B=.86 

r=.60  r=.54 

A=-4.0  A=-5.6 

B=.59  B=  59 

r=.56  r=.64 

A=1.2  A=1.7 

B=.65  B=.75 

r=.34  r=.56 

A=-1.2A=-1.6 

B=.36  B=  65 

500 

r=.65  r=.71 

A=-7.8  A=-14.2 

B=.62  B=.55 

r=.63  r=.51 

A=-ll.l  A=-14. 

B=.59  B=.46 

r=.63  r=.76 

A=2.3  A=2.3 

B=.7  B=.8 

r=.55  i=.59 

A=-1.2  A=-1.6 

B=.58  B=.64 

300 

r=.94  r=.95 

A=-5.8  A=-6.2 

B=.85  B=.81 

r=.59  r=.59 

A=-23.  A=-25. 

B=.52  B=.48 

r=.72  r=.72 

A=3.6  A=5.5 

B=.72  B=.73 

r=.66  r=.67 

A=-2.3  A=-3.4 

B=.71  B.64 

200 

r=.72  r=.81 

A=-21.  A=-20. 

B=.60  B=.81 

r=.73  r=.75 

A=3.8  A=3.5 

B=.77  B=  86 

r=.67  r=.77 

A=3.8  A=3.5 

B=.83  B=.73 

Entire  Data 

r=.97  r=.97 

A=-2.3  A=-1.8 

B=.94  B=.93 

r=.89  r=.87 

A=-2.6  A=-3.9 

B=1.0  B=.94 

r=.70  r=.76 

A=1.7  A=2.4 

B=.77  B=  83 

r=.58  r=.66 

A=-1.8  A=-2.4 

B=.64  B=.67 
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Figure  8A.  Relationship  between  the  GSM 
(top)  for  12-h  forecasting  data  and  upper 
air  sounding  data  for  temperature  for  the 
24-h  forecast  data  (bottom). 


Figure  8C.  Same  as  flgure  8A,  except  for 
east-west  components  of  horizontal  wind 
vectors. 


GSH 


Figure  8D.  Same  as  figure  8A,  except  for 
north-south  components  of  horizontal  wind 
vectors. 


From  the  table  and  figures,  the  following  can  be  inferred: 

•  Among  the  four  parameters,  temperature  shows  the  best  agreement  between  the 
GSM  and  the  sounding  data.  Dew  point  temperature  and  wind  components 
(particularly  the  north-south  wind  component)  show  considerable  dispersion. 

•  For  both  temperature  and  dew  point  temperature,  the  GSM  forecast  data  tends 
to  show  lower  values  than  observed  at  higher  pressure  levels.  See  the  values 
of  A  at  500,  300  and  200-mb  levels.  Restated,  the  GSM  atmosphere  is  drier 
than  the  observed  atmosphere  in  the  upper  levels. 

•  For  the  wind  components,  the  GSM  tends  to  underforecast  wind  speeds  in 
larger  wind  speed  cases  (such  as  jet  streams). 

•  The  GSM  forecasting  skills  for  12  and  24  h  do  not  show  significant  differences. 


•  When  the  GSM  forecast  data  is  used  as  the  initial  and  boundary  condition  data 
for  the  BFM,  it  should  be  improved  as  much  as  possible,  particularly  for  dew 
point  temperature,  horizontal  wind  components. 

The  U.S.  Air  Force  Air  Weather  Service  publication  cites  the  following 

shortcomings  of  the  GSM  forecast,  as  the  findings  of  the  the  U.S.  Air  Force  Global 

Weather  Central:  [34] 

•  Long  waves  are  not  handled  well  because  the  model  was  designed  for  flight- 
level  forecasting. 

•  Forecast  wind  speeds  are  typically  too  slow  below  the  100  mb  level. 

•  Moisture  forecasts  are  typically  too  dry. 


94 
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•  The  model  is  generally  too  warm  in  the  troposphere. 

•  Geopotential  heights  are  generally  too  high. 

From  the  list,  2  and  3  are  in  agreement  with  our  findings;  however,  4  does  not 
agree,  which  shows  the  GSM  temperature  as  slightly  cooler  than  observation  in  the 
lower  levels  and  slightly  warmer  in  higher  levels.  We  cannot  compare  1  and  5 
because  no  study  has  been  done.  Our  experience  has  also  shown  that  the  GSM 
analysis  fields  tend  to  produce  interpolated  grid  point  surface  level  wind  speeds 
that  are  too  high,  and  surface  level  temperatures  that  are  either  too  low  (during 
afternoons)  or  too  high  (during  late  night/early  morning).  A  more  detailed  study 
using  the  data  of  many  different  days  and  areas  is  necessary. 

6.2  Initialization  Fields  for  the  BFM  Forecasting  Calculation 

As  mentioned  previously,  the  BFM  forecast  model  is  initialized  by  nudging  the 
BFM  field  for  3  h  fi*om  Tq  -3  to  Tq  to  the  three-dimensional  fields  produced  by  the 
objective  analysis,  where  Tq  is  the  forecast  base  time.  The  three-dimensional  fields 
are  created  by  compositing  the  GSM  data  and  upper  air  sounding  data 
(composited  data)  when  both  are  available.  Furthermore,  when  surface  observation 
data  is  available,  the  third  and  fourth  layers  of  the  BFM  data  are  modified  by 
nudging  to  the  surface  data  for  3  h  from  To  -  3  to  Tq  In  the  following,  the 
differences  of  the  initialization  data  produced  by  the  above  methods  are  shown, 
using  the  data  over  the  Bosnia  area  at  12  GMT,  24  Aug  96. 

Figure  9  shows  a  horizontal  wind  vector  distribution  at  10  m  AGL,  produced  by 
compositing  the  GSM  00  GMT  base  time  12-h  forecast  field  with  the  available 
upper-air  sounding  data  at  12  GMT,  24  Aug  95.  Wind  directions  are 
predominantly  from  the  west.  Figure  10  shows  observed  horizontal  wind  vectors 
at  the  same  level  and  at  the  same  time.  Observed  wind  directions  at  the  surface 
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were  southerly.  The  BFM  fields  are  nudged  toward  these  fields  for  3  h  from  T  q  - 
3  (9  GMT)  to  To  (12  GMT). 

Figure  11  (top)  is  the  stream  line  distribution  at  the  10  m  level  obtained  by 
nudging  for  3  h  to  both  fields  (the  composited  data  and  the  surface  observation 
data).  For  comparison,  figure  1 1  (bottom)  is  the  streamline  obtained  by  nudging 
to  only  the  composited  data.  Nudging  to  the  surface  observation  data  has  produced 
significantly  different  fields  for  initialization. 

Figure  12  (top)  and  (bottom)  are  the  surface  temperature  distributions  at  the  same 
time.  The  top  figure  is  the  result  by  nudging  to  both  the  composited  data  and  the 
surface  observed  data,  and  (bottom)  is  the  result  obtained  by  nudging  only  to  the 
composited  data.  For  this  case,  there  are  not  significant  differences  in  the  surface 
temperature  distribution  between  the  composited  data  and  the  observation.  The 
distribution  patterns  are  fairly  similar  between  the  two  figures. 

Figure  13  (top)  and  (bottom)  are  the  dew  point  temperature  distributions  at  the 
10-m  level.  The  top  figure  was  produced  by  nudging  to  both  the  composited  data 
and  the  surface  observed  data,  and  the  bottom  figure  by  nudging  to  only  the 
composited  data.  There  are  significant  differences  between  the  two  figures. 
Figure  13  (top)  shows  localized  distributions  of  dew  point  temperature  created  by 
the  strongly  localized  nudging  parameter. 
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Figure  13.  Dew  point  temperature 
distribution  at  the  10-m  level  at  T#  =  0  h 
(top)  for  surface  wind  data  nudged,  and  for 
not  nudged  (bottom). 


6.3  The  Effects  of  Initial  Data  Fields  on  BFM  Output  at  the 
Surface 

The  previous  section  showed  that  the  initialization  fields,  produced  by  nudging  to 
both  the  composited  data  and  the  surface  observation  data,  were  significantly 
different  fi'om  those  produced  by  nudging  to  only  the  composited  data.  To 
examine  how  the  initial  fields  of  the  surface  layers  influence  the  BFM  forecast 
fields,  the  BFM  forecast  calculations  are  made  using  the  two  different  initial  fields. 
Where  surface  nudging  is  used,  the  magnitude  of  the  nudging  is  damped  in  time 
by  a  factor  of  exp(-kt),  where  k  is  an  empirical  coefficient  and  t  is  the  time  after 
To. 

Figure  14  (top  and  bottom)  are  the  streamline  and  isotach  distributions  at  10  m 
AGL  at  Tq  =  3  h.  Figure  14  (top)  is  the  result  obtained  using  both  composited  and 
surface  observation  data,  and  (bottom)  is  the  result  obtained  by  using  only  the 
composited  data.  There  are  significant  differences  between  the  two  images,  which 
show  the  influences  of  the  initial  data  (figure  11).  However,  at  T  q  =  6  h  (figure 
15,  top  and  bottom)  the  two  become  similar  to  each  other.  The  influences  of 
different  initial  surface  data  almost  disappear,  and  the  lateral  boundary  condition 
and  the  model's  boundary  layer  physics  take  control  of  the  surface  layer  wind  flow 
patterns. 

Figures  16  and  17  show  the  surface  temperature  (2  m  AGL)  distributions  at 
forecast  time  Tq  =  3  h  and  T  „  =  6  h,  respectively.  Figures  1 8  and  1 9  are  the  surface 
dew  point  temperature  distributions  at  To  =  3  h  and  Tq  =  6  h,  respectively.  In  these 
figures,  the  top  is  the  result  generated  by  the  composited  and  the  surface 
observation  data,  and  the  bottom  is  the  result  generated  by  the  composited  data 
only.  Similar  to  wind  streamline  distributions,  at  Tq  =  3  h,  the  influences  of  initial 
data  can  be  detected,  but  at  Tq  =  6  h,  the  patterns  of  the  top  and  bottom  images  in 
figures  17  and  19  become  very  similar. 
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It  may  be  concluded  that  after  several  hours  of  forecast  calculation  the  influences 
of  surface  data  disappear  in  the  BFM  output  fields,  and  that  the  model  physics  and 
the  lateral  boundaiy  conditions  become  dominant  to  the  BFM's  fields  of  surface 
wind,  temperature,  and  dew  point  temperature. 
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Figure  18.  Same  as  figure  13,  except  for 
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Figure  19.  Same  as  figure  13,  except  for 


i 


6.4  Examples  of  Fog/Low  Level  Cloud  and  Non-Convective 
Precipitation 

The  BFM  was  applied  to  the  Bosnia  area  for  the  24-h  period  of  12  GMT, 
19  December  through  12  GMT,  20  December  1995.  On  20  December  1995,  the 
air  traffic  over  the  area  was  hampered  by  bad  visibility  and  falling  snow.  The 
model  simulated  cloud  formation  over  the  area  during  the  late  evening  of 

19  December  1995,  in  addition  to  total  cloud  coverage  over  the  entire  model 
domain  from  01  GMT  to  12  GMT,  20  December  1995. 

Figure  20  (top  and  bottom)  show  the  cloud  coverages  at  23  GMT  on  19  December 
1995,  and  00  GMT  on  20  December  1995.  The  figures  show  the  distribution  of 
the  weight  of  liquid  water  over  a  unit  area  (g  /m^).  The  darker  the  color,  the  denser 
the  cloud.  The  figures  may  represent  the  cloud  picture  taken  from  far  above.  As 
can  be  seen,  the  area  covered  by  cloud  increased  in  the  1  h,  and  at  01  GMT  on 

20  December  1995,  the  entire  area  was  covered  by  cloud. 

Figure  21  shows  the  vertical  profiles  of  total  water  content  (g/kg),  cloud  water 
content  (g/kg),  and  relative  humidity  (%),  at  the  center  point  of  the  model  domain. 
The  top  image  represents  00  GMT,  20  December  1995,  and  the  bottom  represents 
01  GMT,  20  December  1995.  At  01  GMT,  the  entire  area  was  covered  by  cloud. 
In  1  h,  the  modeled  total  water  content  increased  significantly,  and  cloud  water 
content  profiles  show  significant  increase  in  the  lower  atmosphere  after  1  h.  The 
relative  humidity  has  also  increased  in  1  h,  and  a  deep  layer  of  the  atmosphere  has 
become  saturated  by  01  GMT.  From  the  figures,  it  can  be  seen  that  a  very  small 
portion  of  total  water  is  converted  into  cloud  water. 

Figure  22  shows  the  distribution  of  precipitation  rate  (mm/h)  at  01  GMT  on 
20  December  1995.  The  precipitation  rate  was  calculated  by  equation  (170).  The 
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figure  shows  the  precipitation  rate  varied  from  0.1  to  1.4  mm/h  over  the  area. 
There  was  no  precipitation  over  the  model  domain  at  00  GMT. 

It  should  be  emphasized  that  fog/cloud  and  precipitation  schemes  in  the  model 
have  yet  to  be  thoroughly  evaluated,  and  the  results  shown  should  be  regarded  as 
preliminary.  Further  study  and  development  will  be  required. 
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Figure  21.  Vertical  profiles  of  total  water 
content  (g/kg),  cloud  water  content  (g/kg),  and 
relative  humidity  (%),  at  the  center  of  model 
domain  (top)  at  00  GMT,  and  (bottom)  at 
01  GMT  20  Dec  95. 


Figure  22.  Precipitation  rate  (mm/h) 
distribution  over  Bosnia,  at  01  GMT, 
20  Dec  95.  Maximum  rate  is  1.4  mm/h,  and 
mimimum  rate  is  .1  mm/h. 


7.  Conclusion 


The  forecasting  period  of  the  BFM  has  been  extended  to  24  h,  and  recently  the 
forecast  schemes  for  fog/cloud  and  as  non-convective  precipitation  have  been 
added.  Other  IMETS  Block  II  software  such  as  the  Atmospheric  Sounding 
Program  (ASP)  and  the  Integrated  Weather  Effects  Decision  Aids  (IWEDA) 
utilize  the  BFM  output  file  to  calculate  the  values  of  their  parameters.  ASP 
generates  products  such  as  probability  of  thunderstorm  occurrence,  icing, 
visibility,  turbulence,  etc. 

The  terrain  elevation  data  production  scheme  is  executed  firom  the  BFM  execution 
GUI  when  the  BFM  is  applied  to  a  new  model  area.  For  a  typical  model  area  of 
500  X  500  km,  meteorological  data  over  an  area  extending  to  1,600  x  1,600  km  is 
used,  so  that  typically  four  DTED  CDROMs  are  required. 

The  three-dimensional  objective  analysis  programs  are  also  executed  firom  the 
BFM  execution  GUI.  Depending  on  the  availability  and  quality  of  the  input  data 
(GSM,  upper  air  sounding,  and  surface  data ),  different  combinations  of  input  data 
are  used  to  initialize  the  BFM.  For  all  input  data  available,  the  BFM  is  initialized 
by  the  composited  data  fields  from  the  GSM  and  upper  air  sounding  data,  plus 
surface  sensor  meteorological  data.  For  a  minimal  requirement,  the  BFM  can  be 
initialized  by  a  single  sounding  profile  within  the  model  domain  and  produce 
forecast  field  at  to+  6  h.  During  the  6-h  forecast  cycle,  the  model  fields  are  nudged 
towards  the  initial  fields  (since  no  GSM  lateral  boundary  conditions  exist)  but  with 
an  exponential  damping  in  time  identical  to  that  used  for  surface  nudging.  The 
BFM  execution  GUI  is  capable  of  selecting  the  combination  of  input  data. 

The  evaluation  of  the  BFM  forecasting  capability  is  ongoing,  and  the  results  will 
be  published  in  the  future.  Forecast  capabilities  of  wind,  temperature,  and  relative 
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humidity,  as  well  as  cloud  coverage  and  non-convective  precipitation  amount 
distribution  will  be  evaluated  by  comparing  the  model  output  and  observed  data. 

The  strengths  of  the  BFM  can  be  summarized  as  follows: 

•  The  BFM  is  numerically  stable  and  fast  in  calculation,  and,  as  long  as 
reasonable  input  data  are  used,  the  BFM  produces  numerically  reasonable 
output.  Currently,  on  the  IMETS  BLOCK  II  computer  (SUN  SPARC  20,  75 
MHz,  256  MB  RAM  single  processor),  the  BFM  requires  about  1.5  to  2  h  of 
computing  time  for  the  24-h  forecast  duration,  using  the  model  configuration 
of  51  X  51  X  16  grid  points  and  10  km  horizontal  grid  spacing. 


•  Detailed  atmospheric  structures  in  the  boundary  layers  can  be  seen.  Sea  and 
land  breeze  circulations,  up-and-down  valley  winds,  and  diurnal  variations  of 
temperature  due  to  surface  heating/cooling  are  well  simulated. 

•  GUI  for  the  model  execution  and  forecast  displays  are  well  designed  and  easy 
to  use.  Graphical  displays  of  output  include  isoline  contours  and  shaded  color 
contour  displays  for  scalar  meteorological  forecast  parameters,  and  streamline 
and  vector  displays  for  wind. 

•  The  BFM  can  be  applied  to  any  area  of  the  world,  as  long  as  the  input  data  for 
initialization  and  lateral  boundary  values  are  available. 

The  presently  known  limitations  of  the  BFM  are  the  following: 

•  Both  hydrostatic  and  Boussinesq  approximations  are  assumed.  Extremely  large 
changes  in  the  large-scale  temperature  field  advected  into  the  BFM  grid  (via 
time-dependent  boundary  values  resulting  from  the  GSM)  over  24  h  could 
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possibly  reduce  the  validity  of  the  Boussinesq  approximation,  based  on  the 
BFM  definition  of  the  basic  state  for  0^  such  as,  the  basic  state  0^  that  is  the 
initially  analyzed  0^  field  generated  from  the  three-dimensional  analysis). 

•  No  cloud  microphysics. 

•  No  convective  cloud  parameterization,  thus  no  ability  to  handle  features 
generated  by  convective  systems,  such  as  meso-fronts  (gust  fronts),  meso 
high/low,  unless  they  are  resolved  by  the  initial  fields. 

•  Rather  crude  specification  of  surface  parameters,  such  as  albedo,  emissivity,  soil 
moisture,  snow,  etc.  These  will  be  immediately  addressed  as  real  or  near-real 
time  data  becomes  readily  available  to  the  IMETS. 

•  Somewhat  time-consuming  to  run  at  resolutions  less  than  5  km. 

•  A  very  coarse  scale  global  model  is  currently  used  to  provide  lateral  boundary 
conditions  to  the  BFM. 

Our  future  plan  for  further  development  includes: 

•  A  BFM  forecast  capability  that  will  be  improved  when  the  output  of  the 
regional  scale  forecast  models,  such  as  the  U.  S.  Air  Force's  Relocatable 
Window  Model  (RWM)  or  the  U.S.  Navy  Operational  Regional  Atmospheric 
Prediction  System  (NORAPS),  is  used  instead  of  the  GSM.  We  will  start 
studying  the  utilization  of  regional  scale  model  oufr)ut  to  initialize  the  BFM  and 
to  provide  its  lateral  boundary  conditions. 
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•  The  inclusion  of  a  cumulus/convective  precipitation  scheme  will  be  made,  by 
using  currently  available  schemes,  such  as  Sun  and  Haines  and  Kuo.  [35],  [36] 
Currently,  we  use  16  vertical  layers  with  the  model  top  7000  m  having  the 
highest  elevation  in  the  model  domain.  To  include  a  convective  precipitation 
scheme,  the  number  of  vertical  layers  and  the  model  depth  will  need  to  be 
increased. 

•  Evaluation  studies  of  cloud/non-convective  precipitation  schemes  by 
qualitatively  comparing  the  model  results  with  satellite  photographs  and 
NEXRAD  estimated  precipitation  patterns,  and  quantitatively  comparing 
precipitation  records. 

•  Coding  the  three-dimensional  analysis  program  of  input  data  more  efficiently, 
to  shorten  the  overall  forecast  computational  time. 
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Acronyms  and  Abbreviations 


AFGWC 

Air  Force  Global  Weather  Center 

AGL 

above  ground  level 

ARL 

Army  Research  Laboratory 

ASL 

above  sea  level 

ASP 

Atmospheric  Sounding  Program 

ATRAD 

atmospheric  radiation 

AWDS 

Automated  Weather  Distribution  System 

BED 

Battlefield  Enviromnent  Directorate 

BFM 

Battlescale  Forecast  Model 

CDROM 

Compact  Disk  Read  Only  Memory 

DMA 

Defense  Mapping  Agency 

DIED 

Digital  Terrain  Elevation  Data 

DTEDl 

Level  1  Digital  Terrain  Elevation  Data 

GMT 

Greenwich  Meridian  Time 

GSM 

Global  Spectral  Model 

GUI 

Graphical  User  Interface 

HOTMAC 

Higher  Order  Turbulence  Model  of  Atmospheric 

Circulation 

IMETS 

Integrated  Meteorological  System 

IWEDA 

Integrated  Weather  Effects  Decision  Aids 

MSL 

mean  sea  level 

NORAPS 

Navy  Operational  Regional  Atmospheric  Prediction 

System 

RWM 

Relocatable  Window  Model 

WSMR 

White  Sands  Missile  Range 
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Appendix 
List  of  Symbols 


tridiagonal  matrix 
expression  by  eq.  (88) 

coefficient  of  a  linear  relationship  y=A  +  Bx 
absorption  of  short-wave  radiation  in  cloud 
parameter  defined  by  eq.  (119) 
parameter  defined  by  eq.  (150) 
column  vector 
Bowen  ratio  ( =  Hj/LE ) 
empirical  constant  ( =  16.6 ) 
expression  by  eq.  (89) 

coefficient  of  a  linear  relationship  y  =  A  +  Bx 
parameter  defined  by  eq.  (134) 
parameter  defined  by  eq.  (133) 
parameter  defined  by  eq.  (120) 
parameter  defined  by  eq.  (151) 
expression  by  eq.  (90) 
nudging  coefficient 
nudging  parameter  for  surface  data 
specific  heat  of  air  at  constant  pressure 

parameter,  I/Cq  is  a  characteristic  time  for  conversion  of  cloud 

droplets  to  precipitation 

parameter  defined  by  eq.  (171) 

empirical  parameter  ( =  0.1) 

constant  ( =  0.01) 

parameter  defined  by  eq.  (121) 


D 

d 

E, 

e 

Ss 

Fb 

Fbf 

Fb^ 

Fo 

F, 
F, 

F2 

F3 

F4 

F„ 

f 

G 


G(z) 

Gfii 

Get 

Gl 


parameter  defined  by  eq.  (127) 

mean  distance 

expression  by  eq.  (83) 

water  vapor  pressure 

saturation  water  vapor  pressure 

parameter  defined  by  eq.  (123) 

temperature  fimetion  defined  by  eq.  (170) 

downward  solar  flux  at  the  top  of  cloud 

parameter  defined  by  eq.  (118) 

expression  by  eq.  (8) 

empirical  constant  (  =  1.8  ) 

right-hand  side  of  eq.  (35) 

empirical  constant  ( =  1.33  ) 

right-hand  side  of  eq.  (36) 

right-hand  side  of  eq.  (46) 

right-hand  side  of  eq.  (44) 

parameter  defined  by  eq.  (169) 

Coriolis  parameter 
ground  heat  flux 

solar  radiation  reaching  the  surface 
probability  density  function  for  cloud  formation 
infrared  flux  profile  in  the  layered  cloud 
downward  long- wave  flux  at  the  cloud  top 
upward  long-wave  flux  at  the  cloud  base 
parameter  defined  by  eq.  (125) 
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Gu 

Go 

Gi 

g 

H 

H 

H(x) 

Hs 

^max 

Jd 

Jmax 

K. 

Ky 

K, 

k 

L 

LE 

Lv 

/ 

m 

n 

P 


parameter  defined  by  eq.  (126) 
parameter  defined  by  eq.  (129) 
parameter  defined  by  eq.  (130) 
acceleration  of  gravity 

material  surface  top  of  the  model  in  z  coordinate 
material  surface  top  of  z*  coordinate 
Heaviside  function 
sensible  heat  flux 

maximum  number  of  grid  points  in  x  direction 
Julian  date 

maximum  number  of  grid  points  in  y  direction 
heat  conductivity  of  soil 

horizontal  eddy  dififusivity  in  x  direction  (  =2c(Ax)(Ay)|  5U/ax|  ) 
horizontal  eddy  difhisivity  in  xy  direction  (  =c(Ax)(Ay){|aU/ax|  + 
|av/ay|}) 

horizontal  eddy  difiusivity  in  y  direction  ( =2c(Ax)(Ay)|  aV/Sy  |  ) 

empirical  parameter  to  determine  the  shape  of  weighting  fimction 

von  Karman  constant 

Obukhov  length 

latent  heat  flux 

latent  heat  of  condensation 

turbulence  length  scale 

fimction  defined  by  eq.  (78) 

time-step 

pressure 
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rate  of  release  of  precipitation 
Pgj  surface  pressure 

turbulence  Prandtl  number 

P,  pressure  at  height 

Pj  pressure  at  height 

^  threshold  value  for  cloud  water  to  produce  precipitation 

Q(c,  cr)F  parameter  defined  by  eq.  ( 1 72) 

Q,  parameter  defined  by  eq.  (160) 

Q,  mixing  ratio  of  liquid  water 

Qj  saturation  mixing  ratio 

Q.  scale  of  moisture 

q  water  vapor  mixing  ratio 

qV2  turbulence  kinetic  energy  of  unit  mass  of  the  air  ( q^=u^  +  v^  +  w^ ) 

R  parameter  defined  by  eq.  ( 1 1 5) 

fraction  of  cloud  coverage 
Rj  gas  constant  of  the  dry  air 

R^  reflection  of  short-wave  in  cloud 

R^  long-wave  radiation  flux 

Re  radius  of  the  earth 

Rif  flux  Richardson  number 

Rig  gradient  Richardson  number 

Ri^  critical  flux  Richardson  number 

Rl,  incoming  long-wave  radiation 

Rl,  outgoing  long- wave  radiation 

Rj  incoming  direct  solar  radiation 
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gas  constant  for  water  vapor 

near-surface  direct  solar  radiation  flux 

incoming  radiation  at  the  top  of  the  atmosphere 

long- wave  incoming  radiation  to  horizontal  surface 

normalized  distance  between  a  grid  point  (i,  j)  and  Nth  GSM  point. 

turbulence  Schmidt  number 

empirical  constant  ( =  0.2 ) 

function  of  flux  Richardson  number,  defined  by  eq.  (49) 

empirical  constant  ( =  0.2 ) 

temperature 

absolute  temperature  at  cloud  top 

absolute  temperature  at  cloud  base 

dew  point  temperature 

groimd  surface  temperature 

soil  surface  temperature 

liquid  water  temperature 

transmission  of  short-wave  radiation  in  cloud 

soil  temperature 

virtual  temperature 

temperature  at  height  ^>i 

temperature  at  height  ^>2 

scale  of  temperature 

dimensionless  longitudinal  distance  defined  by  eq.  (3) 
clock  time 

equation  of  time  defined  by  eq.  (94) 
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Ut 
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u. 


V 

Vg 

V. 


V’ 
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z 

Zb 

Zc 

z* 

Zg 

Zginax2 


'  true  solar  time 

east- west  component  of  horizontal  wind  vector 
dimensionless  latitudinal  distance  defined  by  eq.  (4) 
parameter  defined  by  eq.  (1 1 1) 
east-west  component  of  geostrophic  wind  vector 
east-west  component  of  target  wind  vector 
X  component  of  horizontal  wind  vector  in  the  GSM 
friction  velocity 

north-south  component  of  horizontal  wind  vector 

north-south  component  of  geostrophic  wind 

north-south  component  of  target  wind  vector 

y  component  of  horizontal  wind  vector  in  the  GSM 

vertical  component  of  wind  vector  in  the  Cartesian  coordinate 

total  cloud  liquid  water  content 

vertical  component  of  wind  vector  in  z*  coordinate 

DTED  elevation  data 

average  height 

Cartesian  vertical  coordinate 

cloud  top  height 

cloud  base  height 

vertical  coordinate  used  in  the  BFM 
ground  elevation 

maximum  value  of  the  terrain  elevation  in  the  data  collection 
1600  X  1600  km  domain 
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gmax 


Zo 

Zot 

Zov 

zcrit 


a 


at 
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maximum  value  of  the  terrain  elevation  in  the  BFM  500  x  500  km 
domain 

Cartesian  height  above  sea  level  of  z* 
roughness  length  for  wind 
roughness  length  for  temperature 
roughness  length  for  water  vapor 

model  level  nearest  to  1500  m  AGL  used  to  compute  geostrophic 
winds 

function  of  flux  Richardson  number,  defined  by  eq.  (51) 
angle  of  inclination  to  the  sloped  surface  to  the  horizontal  plane 
coefficient  ( =  0.13) 
coefficient  ( =  0.158) 

backscattered  fraction  of  monodirectional  incident  radiation 

empirical  weight  reduction  factor 

dew  point  depression 

unit  grid  space 

defined  as 

longitudinal  correction  to  time 
declination  of  the  sun 
emissivity  of  the  surface 
non-dimensional  height  ( =  z/L ) 
ratio  of  R/Cp 
potential  temperature 
liquid  potential  temperature 
virtual  potential  temperature 
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fluctuation  of  virtual  potential  temperature 

angle  in  rads  related  to  the  Julian  date  defined  by  eq.  (87) 

longitude  of  the  GSM  grid  point 

smoothing  factor  (0.25  or  0.5) 

wave  length 

parameter  defined  by  eq.  (128) 
parameter  defined  by  eq.  (135) 
parameter  defined  by  eq.  (136) 
function  of  A,  defined  by  eq.  (7) 
zenith  angle 
modified  pressure 
air  density 

Stefan-Boltzman  constant 
parameter  defined  by  eq.  (145) 
parameter  defined  by  eq.  (146) 
parameter  defined  by  eq.  (161) 
optical  thickness  of  cloud 
variable  defined  by  eq.  (114) 
geopotential  height 
latitude 

non-dimensional  wind  profile 
non-dimensional  temperature  profile 
non-dimensional  moisture  profile 
arbitrary  variable 
arbitrary  variable 


empirical  function  to  correct  the  profile  of  temperature 

empirical  function  to  correct  the  profile  of  wind 

azimuth  of  the  projection  of  the  normal  to  the  surface  on  horizontal 

plane 

empirical  function  to  correct  the  profile  of  moisture 
parameter  defined  by  eq.  (16) 
solar  hour  angle 
single  scattering  albedo 
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